Tech Arsenal 1/Tech Arsenal (Arsenal Computer).ISO/tek-03/spbas.zip/SIGPRO.BAS

BASIC Source File | 1993-07-17 | 93KB | 3,544 lines

DECLARE SUB DEGHOST ()
DECLARE SUB HARM (A!(), M(), INV(), S!(), ifset, IFERR)
DECLARE SUB DIMTRANSFORM ()
DECLARE SUB MATCHFILTER ()
DECLARE SUB CROSSCORRELATION ()
DECLARE FUNCTION LN! (X!)
DECLARE SUB SEISMO2 (NEWMODEL!)
DECLARE SUB LEASTDELTA (XGET!(), LEAST!)
DECLARE SUB NEWCONV (ISNP!, IGNP!, S!(), G!(), F!(), CON!(), CC())
DECLARE SUB NEWFOLD (IP!, SIG!(), FOLD!())
DECLARE SUB NEWCONVOLUTION ()
DECLARE SUB SEISMOGRAMMENU ()
DECLARE SUB SEISMOGRAM (N)
DECLARE SUB CONTEST ()
DECLARE SUB FOUR1 (DATQ!(), NN!, ISIGN!)
DECLARE SUB REALFT (DATQ!(), N!, ISIGN!)
DECLARE SUB TWOFFT (DATA1!(), DATA2!(), FFT1!(), FFT2!(), N!)
DECLARE SUB CONVOLVE2 (DATQ!(), N!, RESPNS!(), M!, ISIGN!, ANS!())
DECLARE SUB NUMERICALRECIPESMENU ()
DECLARE SUB CONVOLV (N!, NP2!, S!(), G!())
DECLARE SUB SINCLPF (N!, DELT!, F!())
DECLARE SUB CONVOLUTION ()
DECLARE SUB GETFILTER (FILT!(), FREQ!(), NP2!, FSTITLE$, FTTITLE$)
DECLARE SUB BANDPASSFILTER (S(), TIME(), CN(), FREQ(), RA(), IA(), PHI(), NP2)
DECLARE SUB MENUSIGNALPRO ()
DECLARE SUB GETPOWEROFTWO (NP!, M!, NP2!)
DECLARE SUB RESETALL ()
DECLARE SUB GWHILBERT (HILBERTTRANSFORM!(), W!())
DECLARE SUB FOLDANDSHIFT (IP, XY!(), YX!())
DECLARE SUB SCENTER (OLDCOL!, NEWCOL!, LABLE$)
DECLARE SUB NPRINTFORMAT (NUM!, PRNFMT$)
DECLARE SUB IFEVEN (NUM1!, EVENORODD!)
DECLARE SUB NCENTER (ROW!, COLUMN!, NUMBERLABLE!)
DECLARE SUB WAITFORKEY ()
DECLARE SUB NIFDEFAULT (R$, R!, DEFAULTVAL!)
DECLARE SUB FASTFOURIER (INPUTFILE$, INFILENUM!, OUTPUTFILE$, OUTFILENUM!, NP!)
DECLARE SUB FFT (NP!, M, RA!(), IA!())
DECLARE SUB PAGE ()
DECLARE SUB GRAPHTEST ()
DECLARE SUB GRAPH (SCOUNT!(), S!(), DATATOTAL!, XMIN!, XMAX!, YMIN!, YMAX!, PLOTPOINT%, T$, X$, Y$)
DECLARE SUB GETMINMAX (XGET!(), YGET!(), DATACOUNT!, XMIN!, XMAX!, YMIN!, YMAX!)
DECLARE SUB GETARRAYSIZE (INPUTFILE$, INFILENUM!, XDATATOTAL!, YDATATOTAL!)
DECLARE SUB CNLINEAMPSPEC (A, B, C)
DECLARE SUB XYZDATAFROMARRAY (X, Y, XY())
DECLARE SUB GETINPUTFILE (INPUTFILE$, INFILENUM)
DECLARE SUB GETOUTPUTFILE (OUTPUTFILE$, OUTFILENUM)
DECLARE SUB ANOFCOS (A, B, C, D, E)
DECLARE SUB BNOFSIN (A, B, C, D, E)
DECLARE SUB PHASE (A, B, PHI)
'COMMON SHARED XDATATOTAL, YDATATOTAL, ROWS, COLUMNS
COMMON SHARED INPUTFILE$, OUTPUTFILE$, PAGECOUNT
COMMON SHARED WORKINGDIR$, MAXSUMMATION, MINSUM, FFT2
COMMON SHARED DELTAT, DELTAF, DELTAW, TMAX, FMAX
COMMON SHARED COMMONTITLE$

DEL$ = CHR$(127)
PIE$ = CHR$(227)
PHI$ = CHR$(237)
RHO$ = CHR$(235)

REM $DYNAMIC

REDIM SHARED XYDATA(XDATATOTAL, YDATATOTAL)   'ROWS,COLUMNS
REDIM SHARED XYFOLDED(XDATATOTAL, YDATATOTAL) 'COLUMNS, ROWS
DIM SHARED AMPOFA!(51), AMPOFA2!(51), FOFT!(51), FOFT2!(51)

'SCREEN 12
'WIDTH 80, 60

CONST PI = 3.141593

CONST BLACK% = 0
CONST BLUE% = 1
CONST GREEN% = 2
CONST CYAN% = 3
CONST RED% = 4
CONST MAGENTA% = 5
CONST BROWN% = 6
CONST WHITE% = 7
CONST DARKGREY% = 8
CONST LIGHTBLUE% = 9
CONST LIGHTGREEN% = 10
CONST LIGHTCYAN% = 11
CONST LIGHTRED% = 12
CONST LIGHTMAGENTA% = 13
CONST LIGHTYELLOW% = 14
CONST BRIGHTWHITE% = 15

CONST FALSE = 0, TRUE = NOT FALSE
CLS


RESTART:
DO
COLOR GREEN%
        'DISPLAY MENU
CLS 0

PRINT "WELCOME TO THE WORLD OF GEOPHYSICS PROGRAMS"
PRINT ""
PRINT "PLEASE SELECT YOUR CHOICE FROM THE MENU"
PRINT ""
PRINT "         (1) FOLD AND SHIFT DATA"
PRINT "         (2) NUMERICAL RECIPES"
PRINT "         (3) GRAPH TEST"
PRINT ""
PRINT "         (4) SIGNAL PROCESSING"
PRINT ""
PRINT "         (5) EXIT TO DOS SHELL (TYPE 'EXIT' TO RETURN TO PROGRAM)"
PRINT ""
PRINT "         (6) GET SIGNAL BY TRAPEZOIDAL RULE "
PRINT
PRINT "         (7) "
PRINT ""
PRINT "         (8) NEW CONVOLUTION"
PRINT ""
PRINT "         (9) SEISMOGRAM"
PRINT
PRINT "        (10) QUIT"
INPUT "ENTER A NUMBER TO INDICATE YOUR CHOICE "; CHOICE%
PRINT ""
'RESPOND TO CHOICE
SELECT CASE CHOICE%
        CASE 1:
                CALL GETINPUTFILE(INPUTFILE$, INFILENUM)

                CALL GETOUTPUTFILE(OUTPUTFILE$, OUTFILENUM)

        CASE 2:


        CASE 3: 'GRAPH TEST
                 CALL GRAPHTEST
                 'CALL RESETALL

        CASE 4:

                CALL MENUSIGNALPRO


        CASE 5: COLOR 9
        PRINT "TYPE EXIT TO RETURN TO PROGRAM"
                SHELL
                COLOR 2

        CASE 6:


        CASE 7:

        CASE 8:  CALL NEWCONVOLUTION

        CASE 9:

                CALL SEISMOGRAMMENU


        CASE 10: DONE = TRUE
                END

END SELECT       'SELECT CASE CHOICE
LOOP UNTIL DONE

END

SUB ANOFCOS (N, H, TZERO, PERIOD, A)
'An=(2H/(n*PI) * SIN((2n*PI*TZERO)/PERIOD)  'n=1,2,3...
IF (PERIOD = 0) THEN : A = 0: EXIT SUB
IF N = 0 THEN : A = 4 * H * TZERO / PERIOD: EXIT SUB   'define DC Value=Ao
A = (2 * H / (N * PI)) * SIN((2 * N * PI * TZERO) / PERIOD)

END SUB

REM $STATIC
SUB BANDPASSFILTER (S(), TIME(), CN(), FREQ(), RA(), IA(), PHASEANGLE(), NP2)

REDIM FPOINT(0 TO 3), BPFILTER(0 TO NP2), COUNT(0 TO NP2)

FOR I = 0 TO 3
PRINT "F("; I; ")"; " IS "; : INPUT FPOINT(I)
NEXT I
FOR I = LBOUND(FPOINT) TO UBOUND(FPOINT)
PRINT FPOINT(I)
NEXT I


PRINT "READY TO GRAPH SIGNAL": WAITFORKEY
CALL GETMINMAX(TIME(), S(), NP2, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "INPUT SIGNAL (" + DEL$ + " = " + STR$(DELTAT) + " sec)"
CALL GRAPH(TIME(), S(), NP2, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, "TIME (sec)", "S(t)")

CALL FFT(NP2, M, RA(), IA())  'SEND TIME DOMAIN S()=RA(), GET BACK RA() IN
       ' FREQUENCY DOMAIN, AND IN RECTANGULAR COORDS

FOR N = 0 TO NP2
CALL CNLINEAMPSPEC(RA(N), IA(N), CN(N))  '= |S(w)|,AMPSPEC IN F DOMAIN, NOW POLAR COORDS
CALL PHASE(IA(N), RA(N), PHASEANGLE(N))  '=PHI(w), F DOMAIN, POLAR COORDS
NEXT N

PRINT "READY TO GRAPH AMPLITUDE  SPECTRUM "
CALL GETMINMAX(FREQ(), CN(), NP2, XMIN, XMAX, YMIN, YMAX)
XMAX = FMAX
CALL GRAPH(FREQ(), CN(), NP2, XMIN, XMAX, YMIN, YMAX, 0, "AMPLITUDE SPECTRUM OF THE INPUT SIGNAL", "FREQUENCY (w)", "S(w)")

CALL GWHILBERT(CN(), FREQ())

FOR I = (NP2 * .5) + 1 TO NP2
CN(I) = 0: PHASEANGLE(I) = 0
NEXT I



FOR I = 0 TO NP2 * .5
IF FREQ(I) < FPOINT(0) OR FREQ(I) > FPOINT(3) THEN
BPFILTER(I) = 0
ELSEIF FREQ(I) <= FPOINT(1) THEN
BPFILTER(I) = SIN((2 * PI * (FREQ(I) - FPOINT(0))) / (4 * (FPOINT(1) - FPOINT(0))))
ELSEIF FREQ(I) < FPOINT(2) THEN
BPFILTER(I) = 1
ELSE
BPFILTER(I) = COS((2 * PI * (FREQ(I) - FPOINT(2))) / (4 * (FPOINT(3) - FPOINT(2))))
END IF
CN(I) = CN(I) * BPFILTER(I)
NEXT I


PRINT "READY TO GRAPH THE BAND PASS FILTER": WAITFORKEY
CALL GETMINMAX(FREQ(), BPFILTER(), NP2 * .5, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "BAND PASS FILTER"
XMAX = FMAX
YMIN = 0
CALL GRAPH(FREQ(), BPFILTER(), NP2 * .5, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, "FREQUENCY (Hz)", "BP FILTER")


FOR I = 0 TO NP2    'CONVERT BACK TO RECTANGULAR BEFORE FFT
RA(I) = CN(I) * COS(PHASEANGLE(I))
IA(I) = CN(I) * (SIN(PHASEANGLE(I))) * -1
NEXT I


CALL FFT(NP2, M, RA(), IA())

PRINT "READY FOR THE FILTERED SIGNAL": WAITFORKEY
CALL GETMINMAX(TIME(), RA(), NP2, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "THE BAND PASS FILTERED SIGNAL"
XMAX = TMAX
CALL GRAPH(TIME(), RA(), NP2, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, "TIME (sec)", "S(t)")

END SUB

REM $DYNAMIC
SUB BNOFSIN (N, H, TZERO, PERIOD, B)
IF N = 0 THEN : B = 0: EXIT SUB
B = -H * PERIOD / (N * PI)
B = B * (COS(2 * N * PI * TZERO) - COS(2 * N * PI * (-TZERO)))
END SUB

REM $STATIC
SUB CNLINEAMPSPEC (A, B, C)
C = SQR(A ^ 2 + B ^ 2)
END SUB

SUB CONVOLUTION         'THIS WORKS !!!!!!!!!!!!!!!!!
SNP = 0: GNP = 0: ISNP = 0: IGNP = 0: CONTOTAL = 0: DIVISOR = 1
PRINT "SOLVE FOR h(t) = s(t) * g(t)"
PRINT "FOR THE SIGNAL S(t)"
CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
SINPUTFILE$ = INPUTFILE$: SINFILENUM = INFILENUM
CALL GETOUTPUTFILE(OUTPUTFILE$, OUTFILENUM)
'INPUTFILE$ = "S.DAT"
CALL GETARRAYSIZE(SINPUTFILE$, SINFILENUM, XDATATOTAL, YDATATOTAL)
ISNP = YDATATOTAL - 1: START% = 0

PRINT "FOR THE SIGNAL G(t-T)"
CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
GINPUTFILE$ = INPUTFILE$: GINFILENUM = INFILENUM
CALL GETARRAYSIZE(GINPUTFILE$, GINFILENUM, XDATATOTAL, YDATATOTAL)
IGNP = YDATATOTAL - 1: START% = 0

IF ISNP < IGNP THEN
GNP = IGNP
SNP = GNP
ELSEIF IGNP < ISNP THEN
SNP = ISNP
GNP = SNP
ELSE
SNP = ISNP
GNP = SNP
END IF

PRINT "ISNP "; ISNP; " IGNP "; IGNP; " SNP "; SNP; " GNP "; GNP
INPUT DUM$
REDIM S(0 TO SNP - 1), G(0 TO GNP - 1), GFOLD(0 TO GNP - 1), CON(0 TO 2 * SNP - 1)
REDIM GCOUNT(0 TO GNP - 1), SCOUNT(0 TO SNP - 1)

OPEN SINPUTFILE$ FOR INPUT AS #SINFILENUM
FOR I = 0 TO ISNP - 1
INPUT #SINFILENUM, S(I)
NEXT I
CLOSE #SINFILENUM
OPEN GINPUTFILE$ FOR INPUT AS #GINFILENUM
FOR I = 0 TO IGNP - 1
INPUT #GINFILENUM, G(I)
NEXT I
CLOSE #GINFILENUM

CALL FOLDANDSHIFT(GNP, G(), GFOLD())

IF ISNP < IGNP THEN
FOR I = ISNP + 1 TO GNP - 1
S(I) = 0
NEXT I
ELSEIF IGNP < ISNP THEN
FOR I = IGNP + 1 TO SNP - 1
G(I) = 0
NEXT I
END IF

PRINT "ISNP "; ISNP; " SNP "; SNP
PRINT "IGNP "; IGNP; " GNP "; GNP

FOR I = 0 TO SNP - 1
SCOUNT(I) = I
GCOUNT(I) = I
NEXT I

'*****************************************
LX = SNP: LB = GNP
LY = 0: LY = LX + LB: CONTOTAL = LY
REDIM LX(1 TO LX), LB(1 TO LB), CY(1 TO LY)

PRINT "LX "; LX; " LB "; LB; " LY "; LY
REDIM CONCOUNT(0 TO CONTOTAL)
FOR I = 0 TO CONTOTAL
CONCOUNT(I) = I
NEXT I
FOR I = 1 TO SNP
LX(I) = S(I - 1)
LB(I) = G(I - 1)
NEXT I

FOR I = 1 TO LX
FOR J = 1 TO LB
CY(I + J - 1) = CY(I + J - 1) + LX(I) * LB(J)
NEXT J
NEXT I

FOR I = 0 TO SNP * 2 - 1
CON(I) = CY(I + 1)
NEXT I
ERASE CY, LX, LB

INPUT DUM$
'*****************************************
PRINT "READY TO GRAPH S(t)"
P = SNP
CALL GETMINMAX(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX)

TITLE$ = "S(t)"
YTITLE$ = "S(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH g(t)"
P = GNP
CALL GETMINMAX(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "g(t)"
YTITLE$ = "g(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE FOLDED SIGNAL"
P = GNP
CALL GETMINMAX(GCOUNT(), GFOLD(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "G(t) FOLDED"
YTITLE$ = "G(t) FOLDED"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), GFOLD(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)


PRINT "READY TO GRAPH THE CONVOLUTION"
'LINE INPUT "ENTER THE TITLE "; TITLE$

CALL GETMINMAX(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX)
PRINT "CONTOTAL "; CONTOTAL
PRINT "XMIN "; XMIN
PRINT "XMAX "; XMAX
TITLE$ = "CONVOLUTION p(t)=s(t) * g(t)"
YTITLE$ = "P(t)"
XTITLE$ = ""
TITLE$ = CT$
INPUT DUM$
CALL GRAPH(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)
PRINT "WRITING TO "; OUTPUTFILE$
OPEN OUTPUTFILE$ FOR OUTPUT AS #OUTFILENUM
PRINT #OUTFILENUM, "I", "S(I)", "G(I)", "GFOLD(I)"
FOR I = 0 TO SNP - 1
PRINT #OUTFILENUM, I, ",", S(I), ",", G(I), ",", GFOLD(I)
NEXT I
PRINT #OUTFILENUM, "I", "CON(I)"
FOR I = 0 TO SNP * 2 - 1
PRINT #OUTFILENUM, I, ",", CON(I)
NEXT I
CLOSE #OUTFILENUM

ERASE SCOUNT, GCOUNT, CONCOUNT, S, G, CON, GFOLD



END SUB

SUB CONVOLV (ISNP, NP2, S(), G())
'THIS IS AN OLD VERSION IT DOES NOT WORK ??????????????
IGNP = ISNP: DIVISOR = 1
SNP = NP2: GNP = NP2

REDIM SCOUNT(START% TO NP2), CON(START% TO NP2 * 2)
REDIM GCOUNT(START% TO NP2), GFOLD(START% TO NP2)

FOR I = 0 TO NP2 - 1
SCOUNT(I) = I
GCOUNT(I) = I
NEXT I

CALL FOLDANDSHIFT(NP2, G(), GFOLD())

FOR I = ISNP + 1 TO NP2
GFOLD(I) = 0
S(I) = 0
NEXT I

'********************************
PRINT "STARTED CONVOLUTION AT "; TIME$
SUM = 0
FOR J = 0 TO NP2 - 1
CON(J) = S(J) * GFOLD(NP2 - 1 - J)
SUM = SUM + CON(J)
CON(J) = SUM / DIVISOR
NEXT J
COUNTER = 0
FOR J = NP2 TO (2 * NP2 - 2)
SUM = 0:  COUNTER = COUNTER + 1
FOR I = COUNTER TO NP2 - 1
SUM = S(I) * GFOLD(I - COUNTER) + SUM
NEXT I
CON(J) = SUM / DIVISOR
NEXT J
CON(2 * NP2 - 1) = 0

CONTOTAL = NP2 * 2
REDIM CONCOUNT(0 TO CONTOTAL - 1)
FOR I = 0 TO CONTOTAL - 1
CONCOUNT(I) = I
NEXT I
PRINT "ENDED CONVOLUTION AT "; TIME$
'*********************************
'PRINT "OUTPUT WRITTEN TO "; OUTPUTFILE$
'OPEN OUTPUTFILE$ FOR OUTPUT AS #OUTFILENUM
'OUTFILENUM = FREEFILE
'OPEN "DUM.OUT" FOR OUTPUT AS #OUTFILENUM

'PRINT #OUTFILENUM, "I", "S(I)", "G(I)"
'TMP$ = "############.#######"
'FOR I = 0 TO ISNP - 1
'PRINT #OUTFILENUM, I,
'PRINT #OUTFILENUM, USING TMP$; S(I); G(I)
'NEXT I
'PRINT #OUTFILENUM, "I", "CON(I)"
'FOR I = 0 TO CONTOTAL - 1
'PRINT #OUTFILENUM, I,
'PRINT #OUTFILENUM, USING TMP$; CON(I)
'NEXT I
'CLOSE

PRINT "READY TO GRAPH S(t)"
P = NP2 / 2
CALL GETMINMAX(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "S(t)"
YTITLE$ = "S(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH g(t)"
P = NP2
CALL GETMINMAX(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "g(t)"
YTITLE$ = "g(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE FOLDED SIGNAL"
P = NP2
CALL GETMINMAX(GCOUNT(), GFOLD(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "G(t) FOLDED"
YTITLE$ = "G(t) FOLDED"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), GFOLD(), IGNP, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE CONVOLUTION"
'LINE INPUT "ENTER THE TITLE "; TITLE$
TITLE$ = "Trace1 * f(t), f0 = 40 Hz"
CALL GETMINMAX(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX)
PRINT "CONTOTAL "; CONTOTAL
PRINT "XMIN "; XMIN
PRINT "XMAX "; XMAX
INPUT DUM$
'TITLE$ = "CONVOLUTION p(t)=s(t) * g(t)"
YTITLE$ = "P(t)"
XTITLE$ = ""
CALL GRAPH(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)
'WAITFORKEY
END SUB

SUB CROSSCORRELATION
'NOTE:  SINCE THE CONVOLUTION ROUTINE I'M USING FOLDS THE 2nd SIGNAL, I
'       NEED TO FOLD THE 2nd SIGNAL IN THE CORRELATION ROUTINE BEFORE
'       I SEND IT TO THE CONVOLUTION IN ORDER TO REVERSE THE EFFECTS OF
'       THE FOLDING IN THE CONVOLUTION

SNP = 0: GNP = 0: ISNP = 0: IGNP = 0: CONTOTAL = 0: DIVISOR = 1
PRINT "CROSS CORRELATION"
PRINT "SOLVE FOR h(t) = s(t) * g(t)"
PRINT "FOR THE SIGNAL S(t)"
CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
SINPUTFILE$ = INPUTFILE$: SINFILENUM = INFILENUM
CALL GETOUTPUTFILE(OUTPUTFILE$, OUTFILENUM)
'INPUTFILE$ = "S.DAT"
CALL GETARRAYSIZE(SINPUTFILE$, SINFILENUM, XDATATOTAL, YDATATOTAL)
ISNP = YDATATOTAL - 1
PRINT "ISNP,YT "; ISNP, YDATATOTAL

PRINT "FOR THE SIGNAL G(t-T)"
CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
GINPUTFILE$ = INPUTFILE$: GINFILENUM = INFILENUM
CALL GETARRAYSIZE(GINPUTFILE$, GINFILENUM, XDATATOTAL, YDATATOTAL)
IGNP = YDATATOTAL - 1
PRINT "IGNP, YT "; IGNP, YDATATOTAL

'**************************************

IF ISNP < IGNP THEN
GNP = IGNP
SNP = GNP
ELSEIF IGNP < ISNP THEN
SNP = ISNP
GNP = SNP
ELSE
SNP = ISNP
GNP = SNP
END IF
PRINT "ISNP "; ISNP; " IGNP "; IGNP; " SNP "; SNP; " GNP "; GNP

START% = 0
REDIM S(START% TO SNP), SCOUNT(START% TO SNP), CON(START% TO SNP * 2)
REDIM G(START% TO GNP), GCOUNT(START% TO GNP), GFOLD(START% TO GNP)

FOR I = 0 TO SNP - 1
SCOUNT(I) = I
GCOUNT(I) = I
NEXT I

IF ISNP < IGNP THEN
FOR I = ISNP + 1 TO GNP
S(I) = 0
NEXT I
ELSEIF IGNP < ISNP THEN
FOR I = IGNP + 1 TO SNP
G(I) = 0
NEXT I
END IF

PRINT "ISNP "; ISNP; " SNP "; SNP
PRINT "IGNP "; IGNP; " GNP "; GNP

OPEN SINPUTFILE$ FOR INPUT AS #SINFILENUM
FOR I = 0 TO ISNP - 1
INPUT #SINFILENUM, S(I)    'THE SIGNAL
NEXT I
CLOSE #SINFILENUM

OPEN GINPUTFILE$ FOR INPUT AS #GINFILENUM
FOR I = 0 TO IGNP - 1
INPUT #GINFILENUM, G(I)    'THE SIGNAL
NEXT I
CLOSE #GINFILENUM

COLOR RED%: PRINT "WORKING...": COLOR GREEN%

CALL FOLDANDSHIFT(GNP, G(), GFOLD())

PRINT "I", "S", "G", "GF"
FOR I = 0 TO SNP - 1
PRINT I, S(I), G(I), GFOLD(I)
NEXT I

'*****************************************
LX = SNP: LB = GNP
LY = 0: LY = LX + LB: CONTOTAL = LY
REDIM LX(1 TO LX), LB(1 TO LB), CY(1 TO LY)

PRINT "LX "; LX; " LB "; LB; " LY "; LY
REDIM CONCOUNT(0 TO CONTOTAL)
FOR I = 0 TO CONTOTAL - 1
CONCOUNT(I) = I
NEXT I
FOR I = 1 TO SNP
LX(I) = S(I - 1)
LB(I) = GFOLD(I - 1)
'LB(I) = G(I - 1)
NEXT I

FOR I = 1 TO LX
FOR J = 1 TO LB
CY(I + J - 1) = CY(I + J - 1) + LX(I) * LB(J)
NEXT J
NEXT I
PRINT "IN NEWCONV"
FOR I = 1 TO 2 * SNP
PRINT CY(I)
NEXT I
FOR I = 0 TO SNP * 2 - 1
CON(I) = CY(I + 1)
NEXT I
ERASE CY, LX, LB
'*****************************************

PRINT "ISNP "; ISNP; " SNP "; SNP
PRINT "IGNP "; IGNP; " GNP "; GNP
PRINT "CONTOTAL "; CONTOTAL
INPUT DUM$

PRINT "OUTPUT WRITTEN TO "; OUTPUTFILE$
OPEN OUTPUTFILE$ FOR OUTPUT AS #OUTFILENUM

PRINT #OUTFILENUM, "I", , "S(I)", "G(I)"
TMP$ = "############.#######"
FOR I = 1 TO ISNP
PRINT #OUTFILENUM, I,
PRINT #OUTFILENUM, USING TMP$; S(I); G(I)
NEXT I
PRINT #OUTFILENUM, "I", "CON(I)"
FOR I = 1 TO CONTOTAL
PRINT #OUTFILENUM, I,
PRINT #OUTFILENUM, USING TMP$; CON(I)
NEXT I
CLOSE

PRINT "READY TO GRAPH S(t)"
P = SNP
CALL GETMINMAX(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX)
'INPUT DUM$
TITLE$ = "S(t)"
YTITLE$ = "S(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH g(t)"
P = GNP
CALL GETMINMAX(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "g(t)"
YTITLE$ = "g(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE FOLDED SIGNAL"
P = GNP
CALL GETMINMAX(GCOUNT(), GFOLD(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "G(t) FOLDED"
YTITLE$ = "G(t) FOLDED"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), GFOLD(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)


PRINT "READY TO GRAPH THE CONVOLUTION CROSS CORRELATION"
'LINE INPUT "ENTER THE TITLE "; TITLE$
TITLE$ = "CROSS CORRELATION OF "
INPUT "CROSS CORRELATION OF "; TCC$
TITLE$ = TITLE$ + TCC$

CALL GETMINMAX(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX)
PRINT "CONTOTAL "; CONTOTAL
PRINT "XMIN "; XMIN
PRINT "XMAX "; XMAX
'TITLE$ = "CONVOLUTION p(t)=s(t) * g(t)"
YTITLE$ = "P(t)"
XTITLE$ = ""
CALL GRAPH(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

'WAITFORKEY
ERASE S, SCOUNT, CON, G, GCOUNT
END SUB

SUB DEGHOST
PRINT "FOR THE SIGNAL E(t)"
CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
SINPUTFILE$ = INPUTFILE$: SINFILENUM = INFILENUM
CALL GETOUTPUTFILE(OUTPUTFILE$, OUTFILENUM)
'INPUTFILE$ = "S.DAT"
CALL GETARRAYSIZE(SINPUTFILE$, SINFILENUM, XDATATOTAL, YDATATOTAL)
ISNP = YDATATOTAL - 1: START% = 0
CALL GETPOWEROFTWO(ISNP, M, NP2)

REDIM S(0 TO ISNP - 1), G(0 TO ISNP - 1), GFOLD(0 TO ISNP - 1), CON(0 TO 2 * NP2 - 1)
REDIM GCOUNT(0 TO ISNP - 1), SCOUNT(0 TO ISNP - 1)
REDIM RA(0 TO NP2 - 1), IA(0 TO NP2 - 1), OMEGA(0 TO NP2 - 1), FREQ(0 TO NP2 - 1)
REDIM TIME(0 TO NP2 - 1), PHASEANGLE(0 TO NP2 - 1), F(0 TO NP2 - 1)
REDIM FP(0 TO NP2 - 1), CN(0 TO NP2 - 1)


OPEN SINPUTFILE$ FOR INPUT AS #SINFILENUM
FOR I = 0 TO ISNP - 1
INPUT #SINFILENUM, S(I)
NEXT I
CLOSE #SINFILENUM


INPUT "ENTER DEPTH OF SOURCE       "; DEPTHOFSOURCE
INPUT "ENTER THE VELOCITY OF WATER "; VH20
INPUT "ENTER Ro                    "; RO

TAUZERO = 2 * DEPTHTOSOURCE / VH20
DELTAT = .01
DELTAF = 1 / (NP2 * DELTAT)


FOR I = 0 TO NP2 - 1
TIME(I) = I * DELTAT
FREQ(I) = I * DELTAF
OMEGA(I) = FREQ(I) * 2 * PI
RA(I) = (1 + RO * COS(TAUZERO * OMEGA(I))) ^ 2
PRINT "RA "; RA(I)  ': INPUT DUM$
IA(I) = RO * SIN(TAUZERO * OMEGA(I))
CALL CNLINEAMPSPEC(RA(I), IA(I), CN(I))  '= |S(w)|,AMPSPEC IN F DOMAIN, NOW POLAR COORDS
CALL PHASE(IA(I), RA(I), PHASEANGLE(I))  '=PHI(w), F DOMAIN, POLAR COORDS
'PHASEANGLE(I) = PHASEANGLE(I) + OMEGA(I) * TAUZERO
PRINT "RA "; RA(I)
PRINT "IA "; IA(I)
PRINT "PHI "; PHASEANGLE(I)  ': INPUT DUM$
F(I) = 1 / (1 + RA(I) ^ 2 + IA(I) ^ 2) ^ .5 * PHASEANGLE(I)
FP(I) = 0
NEXT I


PRINT "READY TO GRAPH THE F"
CALL GETMINMAX(FREQ(), F(), NP2, XMIN, XMAX, YMIN, YMAX)
'XMAX = FMAX
'YMIN = 0
INPUT "XMIN "; XMIN
INPUT "XMAX "; XMAX
INPUT DUM$
CALL GRAPH(FREQ(), F(), NP2, XMIN, XMAX, YMIN, YMAX, 0, FTTITLE$, "FREQUENCY (Hz)", "FILTER")


CALL FFT(NP2, M, F(), FP())    'SENDING F DOMAIN, WILL RETURN T DOMAIN


PRINT "READY TO GRAPH THE F IN TIME"
CALL GETMINMAX(TIME(), F(), NP2, XMIN, XMAX, YMIN, YMAX)
'XMAX = FMAX
'YMIN = 0
INPUT "XMIN "; XMIN
INPUT "XMAX "; XMAX

INPUT DUM$
CALL GRAPH(TIME(), F(), NP2, XMIN, XMAX, YMIN, YMAX, 0, FTTITLE$, "FREQUENCY (Hz)", "FILTER")




END SUB

SUB DIMTRANSFORM
REM $DYNAMIC
'DEFINT I-N
FORWARD = 1
REVERSE = -1

'********************* 2 D TRANSFORM
REM $DYNAMIC
'DEFINT I-N
FORWARD = 1
REVERSE = -1

PRINT "FOR THE FILE CONTAINING THE MULTIDIMENSIONAL DATA "

CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
SINPUTFILE$ = INPUTFILE$: SINFILENUM = INFILENUM
CALL GETOUTPUTFILE(OUTPUTFILE$, OUTFILENUM)
INPUTFILE$ = "S.DAT"
CALL GETARRAYSIZE(SINPUTFILE$, SINFILENUM, XDATATOTAL, YDATATOTAL)
ISNP = YDATATOTAL - 1
PRINT "ISNP,YT "; ISNP, YDATATOTAL
ROW = YDATATOTAL: COL = XDATATOTAL



NDIM = XDATATOTAL - 1
XCOL = XDATATOTAL
YROW = YDATATOTAL
SNDAT = NDIM * YDATATOTAL * 2

'CALL GETPOWEROFTWO(SNDAT, Z, SNP2)
SNP2 = SNDAT

REDIM NN(3)
PROD = 1
FOR I = 1 TO 3
NN(I) = 2 * (2 ^ I)
PROD = PROD * NN(I)
NEXT I



PRINT " PROD IS "; PROD
PRINT "X COL = "; COL; "Y  ROWS = "; ROW
PRINT "SNDAT "; SNDAT; " SNP2 "; SNP2: INPUT DUM$

REDIM DAT(0 TO XDATATOTAL - 1, 0 TO YDATATOTAL - 1), SS(SNP2), INV(SNP2)
'REDIM AMP(0 TO (SNP2 - 1) * PROD), AMPCOUNT(0 TO (SNP2 - 1) * PROD)
'REDIM NEWAMP(0 TO (SNP2 - 1) * PROD)

AMAX = 32767 - 17000           ''  ARRAY MAX 32767
REDIM AMP(0 TO AMAX), AMPCOUNT(0 TO AMAX), NEWAMP(0 TO AMAX)



INF = FREEFILE
OPEN SINPUTFILE$ FOR INPUT AS #INF
AMPCOUNT = 0
FOR ROW = 0 TO YROW - 1
FOR COL = 0 TO XCOL - 2
INPUT #INF, DAT(COL, ROW)
AMPCOUNT(AMPCOUNT) = AMPCOUNT
AMP(AMPCOUNT) = DAT(COL, ROW)
AMPCOUNT = AMPCOUNT + 1
NEXT COL
NEXT ROW

'CALL GETMINMAX(AMPCOUNT(), AMP(), AMPCOUNT, XMIN, XMAX, YMIN, YMAX)
'CALL GRAPH(AMPCOUNT(), AMP(), AMPCOUNT, XMIN, XMAX, YMIN, YMAX, 0, T$, X$, Y$)


'CALL HARM (a(), M%(), INV%(), S(), ifset%, IFERR%)
ISIGN = FORWARD
CALL HARM(AMP(), NN(), INV(), SS(), ISIGN, IFFERR)

NEWAMPCOUNT = 0
FOR I = 0 TO AMPCOUNT - 2
NEWAMP(I) = SQR(AMP(I) ^ 2 + AMP(I + 1) ^ 2)
NEXT I


CALL GETMINMAX(AMPCOUNT(), NEWAMP(), AMPCOUNT, XMIN, XMAX, YMIN, YMAX)
CALL GRAPH(AMPCOUNT(), NEWAMP(), AMPCOUNT, XMIN, XMAX, YMIN, YMAX, 0, T$, X$, Y$)

ERASE NN, DAT, INV, SS, AMP



'**************************************************
END SUB

SUB FASTFOURIER (INPUTFILE$, INFILENUM, OUTPUTFILE$, OUTFILENUM, NP)
CALL GETPOWEROFTWO(NP, M, NP2)
INPUT "ENTER THE SAMPLE RATE IN MSEC (1.0) "; R$: CALL NIFDEFAULT(R$, DELTAT, 1)

INPUT "DO YOU WANT TO ENTER A PHASE SHIFT "; YORN$
ISPHASESHIFT = FALSE
IF UCASE$(YORN$) = "Y" THEN
ISPHASESHIFT = TRUE
INPUT "ENTER THE SHIFT ( + T OR - T) IN msec "; PSHIFT
PSHIFT = PSHIFT * .001   'CONVERT MSEC TO SEC
END IF
INPUT "DO YOU WANT TO FILTER "; YORN$              'FILTER
DOFILT = FALSE
IF UCASE$(YORN$) = "Y" THEN
DOFILT = TRUE
REDIM FILT(0 TO NP2)
END IF

DELTAT = DELTAT * .001      'CONVERT MSEC TO SEC
PRINT "DELTAT = "; DELTAT
PRINT "TOTAL OF S() "; NP
PRINT "NP2 IS "; NP2
PRINT "M = "; M
DELTAF = 1 / NP2 / DELTAT
DELTAW = 2 * PI * DELTAF
PHASESHIFT$ = STR$(PSHIFT)
FMAX = 1 / (2 * DELTAT)
TMAX = NP * DELTAT

PRINT "WORKING...": START% = 0
REDIM AN(START% TO NP2), BN(START% TO NP2)
REDIM CN(START% TO NP2), TIME(START% TO NP2)
REDIM PHASEANGLE(START% TO NP2), SIGNEW(START% TO NP2)
REDIM S(START% TO NP2), FREQ(START% TO NP2), ENV(START% TO NP2)
REDIM RA(START% TO NP2), IA(START% TO NP2), COUNT(START% TO NP2)
IF ISPHASESHIFT = TRUE THEN
REDIM OMEGA(START% TO NP2)
END IF

OPEN INPUTFILE$ FOR INPUT AS #INFILENUM
FOR I = 0 TO NP - 1
INPUT #INFILENUM, S(I)    'THE SIGNAL
RA(I) = S(I)              'SENDING TIME DOMAIN
NEXT I
CLOSE #INFILENUM
FOR I = 0 TO NP2 - 1
TIME(I) = I * DELTAT
FREQ(I) = I * DELTAF
COUNT(I) = I
NEXT I

PRINT "READY TO GRAPH SIGNAL"
CALL GETMINMAX(COUNT(), S(), NP2, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "INPUT SIGNAL (" + DEL$ + " = " + STR$(DELTAT) + " sec)"
CALL GRAPH(COUNT(), S(), NP2, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, "TIME (sec)", "S(t)")

CALL FFT(NP2, M, RA(), IA())  'SEND TIME DOMAIN S()=RA(), GET BACK RA() IN
       ' FREQUENCY DOMAIN, AND IN RECTANGULAR COORDS

IF NOT ISPHASESHIFT THEN   '<> TRUE THEN
FOR N = 0 TO NP2
CALL CNLINEAMPSPEC(RA(N), IA(N), CN(N))  '= |S(w)|,AMPSPEC IN F DOMAIN, NOW POLAR COORDS
CALL PHASE(IA(N), RA(N), PHASEANGLE(N))  '=PHI(w), F DOMAIN, POLAR COORDS
NEXT N
'PRINT "NO PHASESHIFT ": WAITFORKEY
ELSE
FOR N = 0 TO NP2
CALL CNLINEAMPSPEC(RA(N), IA(N), CN(N))  '= |S(w)|,AMPSPEC IN F DOMAIN, NOW POLAR COORDS
CALL PHASE(IA(N), RA(N), PHASEANGLE(N))  '=PHI(w), F DOMAIN, POLAR COORDS
OMEGA(N) = FREQ(N) * 2 * PI
PHASEANGLE(N) = PHASEANGLE(N) + OMEGA(N) * PSHIFT
NEXT N
'PRINT "YES PHASESHIFT ": WAITFORKEY
END IF

FOR I = (NP2 * .5) + 1 TO NP2
CN(I) = 0: PHASEANGLE(I) = 0
NEXT I

PRINT "READY TO GRAPH AMPLITUDE  SPECTRUM "
CALL GETMINMAX(FREQ(), CN(), NP2, XMIN, XMAX, YMIN, YMAX)
XMAX = FMAX
CALL GRAPH(FREQ(), CN(), NP2, XMIN, XMAX, YMIN, YMAX, 0, "AMPLITUDE SPECTRUM OF THE INPUT SIGNAL", "FREQUENCY (w)", "S(w)")

CALL GWHILBERT(CN(), FREQ())

'*********************  THE FILTER
IF DOFILT THEN
'FPOINT(0) = 60: FPOINT(1) = 90: FPOINT(2) = 120: FPOINT(3) = 175

CALL GETFILTER(FILT(), FREQ(), NP2, FSTITLE$, FTTITLE$)

FOR I = 0 TO NP2
CN(I) = CN(I) * FILT(I)
NEXT I


PRINT "READY TO GRAPH THE FILTER"
CALL GETMINMAX(FREQ(), FILT(), NP2 * .5, XMIN, XMAX, YMIN, YMAX)
XMAX = FMAX
YMIN = 0
CALL GRAPH(FREQ(), FILT(), NP2 * .5, XMIN, XMAX, YMIN, YMAX, 0, FTTITLE$, "FREQUENCY (Hz)", "FILTER")

'*****************
REDIM FQ(0 TO NP2), F2(0 TO NP2), F2R(0 TO NP2), F2I(0 TO NP2)
REDIM FPHASE(0 TO NP2)
FOR I = 0 TO NP2
FQ(I) = FREQ(I)
F2(I) = FILT(I)
NEXT I
CALL GWHILBERT(F2(), FQ())
FOR I = 0 TO NP2
F2R(I) = F2(I) * COS(FPHASE(I))
F2I(I) = F2(I) * SIN(FPHASE(I)) * -1
NEXT I

CALL FFT(NP2, M, F2R(), F2I())
DUM = FREEFILE
OPEN "FDUM.OUT" FOR OUTPUT AS #DUM
FOR I = 0 TO NP2
PRINT #DUM, F2(I)
NEXT I
CLOSE #DUM
PRINT "READY TO GRAPH THE FILTER IN TIME DOMAIN"
CALL GETMINMAX(TIME(), F2R(), NP2, XMIN, XMAX, YMIN, YMAX)
PRINT "TMAX "; TMAX
INPUT "ENTER XMAX "; XMAX
'XMAX = NP2 / 2
CALL GRAPH(TIME(), F2R(), NP2, XMIN, XMAX, YMIN, YMAX, 0, FTTITLE$, "TIME", "FILTER")
INPUT DUM$
ERASE FQ, F2, F2R, F2I, FPHASE
END IF
'************************


FOR I = 0 TO NP2    'CONVERT BACK TO RECTANGULAR BEFORE FFT
RA(I) = CN(I) * COS(PHASEANGLE(I))
IA(I) = CN(I) * (SIN(PHASEANGLE(I))) * -1
NEXT I

CALL FFT(NP2, M, RA(), IA())    'SENDING F DOMAIN, WILL RETURN T DOMAIN
FOR I = 0 TO NP2
ENV(I) = (IA(I) ^ 2 + RA(I) ^ 2) ^ .5
NEXT I

'PRINT "READY TO GRAPH HILBERT TRANSFORM"
'HILBERT TRANSFORM = IA(I)   'THE SIGNAL = RA(I)
'CALL GETMINMAX(TIME(), IA(), NP2, XMIN, XMAX, YMIN, YMAX)
'XMAX = TMAX
'IF DOFILT THEN
'TITLE$ = "HILBERT TRANSFORM OF FILTERED SIGNAL (IA)" + FILTERPOINTS$
'ELSE
'TITLE$ = "HILBERT TRANSFORM (IA)"
'END IF
'CALL GRAPH(TIME(), IA(), NP2, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, "TIME (SEC)", "HT(t)",TSCOUNT,)

CALL GETMINMAX(COUNT(), RA(), NP2, XMIN, XMAX, YMIN, YMAX)

IF DOFILT THEN

ELSE
FSTITLE$ = "THE NEW SIGNAL"
END IF

CALL GRAPH(COUNT(), RA(), NP2, XMIN, XMAX, YMIN, YMAX, 0, FSTITLE$, "TIME (SEC)", "S(t)")

'PRINT "READ TO GRAPH ENVELOPE"
'CALL GETMINMAX(TIME(), ENV(), NP2, XMIN, XMAX, YMIN, YMAX)
'XMAX = TMAX
'CALL GRAPH(TIME(), ENV(), NP2, XMIN, XMAX, YMIN, YMAX, 0, "ENVELOPE", "TIME (SEC)", "ENV(t)",TSCOUNT,)

COLOR LIGHTRED%: PRINT "WRITING OUTPUT TO "; OUTPUTFILE$: COLOR GREEN%
OPEN OUTPUTFILE$ FOR OUTPUT AS #OUTFILENUM
PRINT #OUTFILENUM, "N", "TIME(N)", "S(N)", "NEW SIGNAL", "FREQ(N)", "Cn", "HILBERT", "ENVELOPE"
FOR N = 0 TO NP2
CALL NPRINTFORMAT(N, TMP$): PRINT #OUTFILENUM, USING TMP$; N;
PRINT #OUTFILENUM, ",";
CALL NPRINTFORMAT(TIME(N), TMP$): PRINT #OUTFILENUM, USING TMP$; TIME(N);
PRINT #OUTFILENUM, ",";
CALL NPRINTFORMAT(S(N), TMP$): PRINT #OUTFILENUM, USING TMP$; S(N);
PRINT #OUTFILENUM, ",";
CALL NPRINTFORMAT(RA(N), TMP$): PRINT #OUTFILENUM, USING TMP$; RA(N);
PRINT #OUTFILENUM, ",";
CALL NPRINTFORMAT(FREQ(N), TMP$): PRINT #OUTFILENUM, USING TMP$; FREQ(N);
PRINT #OUTFILENUM, ",";
CALL NPRINTFORMAT(CN(N), TMP$): PRINT #OUTFILENUM, USING TMP$; CN(N);
PRINT #OUTFILENUM, ",";
CALL NPRINTFORMAT(IA(N), TMP$): PRINT #OUTFILENUM, USING TMP$; IA(N);
PRINT #OUTFILENUM, ",";
CALL NPRINTFORMAT(ENV(N), TMP$): PRINT #OUTFILENUM, USING TMP$; ENV(N)
NEXT N
CLOSE #OUTFILENUM

END SUB

REM $STATIC
SUB FFT (NP2, M, RA(), IA())

PRINT "BEGINNING FFT AT    "; TIME$

     REM  FFT AFTER COOLEY AND TUKEY
     NV = NP2 / 2
     NM = NP2 - 1
     J = 1
     FOR I = 1 TO NM
     IF I >= J GOTO 2100
     IT = IA(J)
     RT = RA(J)
     RA(J) = RA(I)
     IA(J) = IA(I)
     RA(I) = RT
     IA(I) = IT
2100 K = NV
2110 IF K >= J GOTO 2150
     J = J - K
     K = K / 2
     GOTO 2110
2150 J = J + K
     NEXT I
     FOR L = 1 TO M
     LE = 2 ^ L
     L1 = LE / 2
     RU = 1!
     IU = 0!
     RW = COS(PI / L1)
     IW = SIN(PI / L1)
     FOR J = 1 TO L1
     FOR I = J TO NP2 STEP LE
     IP = I + L1
     RT = RA(IP) * RU - IA(IP) * IU
     IT = IA(IP) * RU + RA(IP) * IU
     RA(IP) = RA(I) - RT
     IA(IP) = IA(I) - IT
     RA(I) = RA(I) + RT
     IA(I) = IA(I) + IT
     NEXT I
     Q1 = IU * IW
     Q2 = RU * RW
     Q3 = IU * RW
     Q4 = RU * IW
     RU = Q2 - Q1
     IU = Q3 + Q4
     NEXT J
     NEXT L

PRINT "ENDED FFT AT    "; TIME$
END SUB

SUB FOLDANDSHIFT (IP, SIG(), FOLD())
'FOLD THE SIGNAL
FOR I = 0 TO IP - 1
FOLD(I) = SIG(IP - 1 - I)
NEXT I
END SUB

SUB GETARRAYSIZE (INPUTFILE$, INFILENUM, XDATATOTAL, YDATATOTAL)
'GET NUMBER OF COLUMNS (TO GET XDATA)
XDATATOTAL = 0: YDATATOTAL = 0: COMMAS = 0: ROWLENGTH = 0
ROWTOTAL = 0: COLUMNTOTAL = 0: ROWS = 0: COLUMNS = 0
OPEN INPUTFILE$ FOR INPUT AS #INFILENUM
LINE INPUT #INFILENUM, ROWLENTH$
ROWLENGTH = LEN(ROWLENTH$)

FOR I = 1 TO ROWLENGTH
IF (MID$(ROWLENTH$, I, 1) = ",") THEN
COMMAS = COMMAS + 1
END IF
NEXT I

XDATATOTAL = COMMAS + 1
COLUMNTOTAL = XDATATOTAL

CLOSE #INFILENUM

'GET NUMBER OF ROWS (TO GET YDATA)
OPEN INPUTFILE$ FOR INPUT AS #INFILENUM

DO
IF EOF(INFILENUM) THEN :  EXIT DO
LINE INPUT #INFILENUM, DUM$
YDATATOTAL = YDATATOTAL + 1
LOOP UNTIL EOF(INFILENUM)

ROWTOTAL = YDATATOTAL
CLOSE #INFILENUM

'INPUT DATA INTO ARRAY

OPEN INPUTFILE$ FOR INPUT AS #INFILENUM
REDIM ROWCOLUMN(ROWTOTAL, COLUMNTOTAL)   'Z=(ROW,COLUMN) = Z=(YDATA,XDATA)
REDIM XYDATA(COLUMNTOTAL, ROWTOTAL)

FOR ROW = 1 TO ROWTOTAL   'YDATA
FOR COLUMN = 1 TO COLUMNTOTAL   'XDATA
IF NOT EOF(INFILENUM) THEN
INPUT #INFILENUM, ROWCOLUMN(ROW, COLUMN)  'ROWCOLUMN(YDATA,XDATA)
END IF
NEXT COLUMN  'XDATA
NEXT ROW     'YDATA

ROWS = XDATATOTAL
COLUMNS = YDATATOTAL
CLOSE #INFILENUM

PRINT "XDATATOTAL "; XDATATOTAL
PRINT "YDATATOTAL "; YDATATOTAL


END SUB

SUB GETFILTER (FILT(), FREQ(), NP2, FSTITLE$, FTTITLE$)
DO
PRINT "SELECT THE FILTER TYPE"
PRINT : PRINT "1) LOW PASS"
PRINT : PRINT "2) HIGH PASS"
PRINT : PRINT "3) BAND PASS"
PRINT : INPUT CHOICE%
LOOP WHILE CHOICE% < 1 OR CHOICE% > 3

SELECT CASE CHOICE%
'FILT(I) = SIN((2 * PI * (FREQ(I) - F0)) / (4 * (F1 - F0)))
'FILT(I) = COS((2 * PI * (FREQ(I) - F2)) / (4 * (F3 - F2)))

CASE 1: ' LOW PASS
INPUT "F2 IS "; F2
INPUT "F3 IS "; F3

FOR I = 0 TO NP2
IF FREQ(I) < F2 THEN
FILT(I) = 1
ELSEIF FREQ(I) <= F3 THEN
FILT(I) = COS((2 * PI * (FREQ(I) - F2) / (4 * (F3) - F2)))
ELSE FILT(I) = 0
END IF
NEXT I

FSTITLE$ = "LOW PASS FILTERED SIGNAL" + STR$(F2) + " " + STR$(F3)
FTTITLE$ = "LOW PASS FILTER" + STR$(F2) + " " + STR$(F3)

CASE 2:  'HIGH PASS
INPUT "F0 IS "; F0
INPUT "F1 IS "; F1

FOR I = 0 TO NP2
IF FREQ(I) < F0 THEN
FILT(I) = 0
ELSEIF FREQ(I) < F1 THEN
FILT(I) = SIN((2 * PI * (FREQ(I) - F0)) / (4 * (F1 - F0)))
ELSE
FILT(I) = 1
END IF
NEXT I

FSTITLE$ = "HIGH PASS FILTERED SIGNAL" + STR$(F0) + " " + STR$(F1)
FTTITLE$ = "HIGH PASS FILTER" + STR$(F0) + " " + STR$(F1)

CASE 3:
'INPUT "ENTER  (1.0) "; R$: CALL NIFDEFAULT(R$, DELTAT, 1)

INPUT "F0 IS (60)  "; R$: CALL NIFDEFAULT(R$, F0, 60)
INPUT "F1 IS (90)  "; R$: CALL NIFDEFAULT(R$, F1, 90)
INPUT "F2 IS (120) "; R$: CALL NIFDEFAULT(R$, F2, 120)
INPUT "F3 IS (175) "; R$: CALL NIFDEFAULT(R$, F3, 175)

FOR I = 0 TO NP2 * .5
IF FREQ(I) < F0 OR FREQ(I) > F3 THEN
FILT(I) = 0
ELSEIF FREQ(I) < F1 THEN
FILT(I) = SIN((2 * PI * (FREQ(I) - F0)) / (4 * (F1 - F0)))
ELSEIF FREQ(I) <= F2 THEN
FILT(I) = 1
ELSE
FILT(I) = COS((2 * PI * (FREQ(I) - F2)) / (4 * (F3 - F2)))
END IF
'CN(I) = CN(I) * FILT(I)
NEXT I

FPOINTS$ = STR$(F0) + " " + STR$(F1) + " " + STR$(F2) + " " + STR$(F3)
FTTITLE$ = "BAND PASS FILTER" + FPOINTS$
FSTITLE$ = "BAND PASS FILTERED SIGNAL " + FPOINTS$

END SELECT
END SUB

SUB GETINPUTFILE (INPUTFILE$, INFILENUM)

'''INPUT "ENTER WORKING DIRECTORY "; WORKINGDIR$
''FILES WORKINGDIR$ + "\*.*"
FILES "*.DAT"
PRINT "ENTER INPUT FILE NAME AND PATH IF NOT IN DEFAULT DIRECTORY "
INPUT INPUTFILE$
INPUTFILE$ = WORKINGDIR$ + INPUTFILE$
PRINT
INFILENUM = FREEFILE
OPEN INPUTFILE$ FOR INPUT AS #INFILENUM: CLOSE #INFILENUM
END SUB

SUB GETMINMAX (XGET(), YGET(), DATACOUNT, XMIN, XMAX, YMIN, YMAX) STATIC
XMIN = 0: XMAX = 0: YMIN = 0: YMAX = 0

PRINT "BEGAN GET GRAPH PARAMS AT    "; TIME$

FOR I = (LBOUND(XGET)) TO (UBOUND(XGET))     'GET MIN'S AND MAX'S
IF XGET(I) < XMIN THEN XMIN = XGET(I)
IF XGET(I) > XMAX THEN XMAX = XGET(I)
NEXT I
FOR I = (LBOUND(YGET)) TO (UBOUND(YGET))
IF YGET(I) < YMIN THEN YMIN = YGET(I)
IF YGET(I) > YMAX THEN YMAX = YGET(I)
NEXT I


PRINT "XMIN "; XMIN
PRINT "XMAX "; XMAX
PRINT "YMIN "; YMIN
PRINT "YMAX "; YMAX
'WAITFORKEY
PRINT "ENDED GET GRAPH PARAMS AT    "; TIME$

END SUB

SUB GETOUTPUTFILE (OUTPUTFILE$, OUTFILENUM)

FILES

PRINT "ENTER OUTPUT FILE NAME AND PATH IF NOT IN DEFAULT DIRECTORY "
INPUT OUTPUTFILE$
OUTPUTFILE$ = WORKINGDIR$ + OUTPUTFILE$
PRINT
OUTFILENUM = FREEFILE

OPEN OUTPUTFILE$ FOR OUTPUT AS #OUTFILENUM: CLOSE #OUTFILENUM
END SUB

SUB GETPOWEROFTWO (NP, M, NP2)
'SUB GETPOWEROFTWO (NP, M, NP2)
M = 0
COUNT = 2
DO
ITESTVAL = 2 ^ COUNT
IF (ITESTVAL >= NP) AND (NP >= 4) THEN
        NP2 = ITESTVAL
        'M = COUNT
        M = COUNT + 1
        NP2 = 2 ^ M
        EXIT DO
END IF
COUNT = COUNT + 1
LOOP WHILE ITESTVAL <= NP
END SUB

SUB GRAPH (X(), Y(), DATATOTAL, XMIN, XMAX, YMIN, YMAX, PLOTPOINT%, TITLE$, XAXIS$, YAXIS$) STATIC
'PLOTPOINT% :   0 = NO POINTS, > 0 YES POINTS , > 1 SPIKE
SCREEN 12: WIDTH 80, 60: CLS 0:  XF = 1: YF = 4 / 3  'SF = 4 / 3:
FORECOLOR% = CYAN%: COLOR FORECOLOR%: XSTOP = 0
CIRCLERADIUS = (.0045 * (XMAX - XMIN))
LXTOPX% = 0: LYTOPY% = 1: PXTOLX% = 2: PYTOLY% = 3
'X AXIS FROM COL 11 TO COL 73 FOR A TOTAL OF 63 COLUMNS (504 PIX)
XEND = XMAX: YEND = YMAX

NUMBEROFXTICKLABLES = 63 / 9   '= 7 X DATA LABLES, EACH 9 COLUMNS WIDE
YINTERVAL% = 9             '45 ROWS (9 Y LABLES)
'XINTERVAL = NUMBEROFXTICKLABLES
XSTEP = (XMAX - XMIN) / (NUMBEROFXTICKLABLES - 1) '=LOGICAL PIXELS PER COLUMN
YSTEP = (YMAX - YMIN) / YINTERVAL%

LXPIXELTOCOL = 63 / (XMAX - XMIN)            '63=504/8
X1 = XMIN: Y1 = YMAX           'UPPER LEFT
X2 = XMAX: Y2 = YMIN           'LOWER RIGHT
Y1 = Y1 * YF: Y2 = Y2 * YF

' 1 TO 58 ROWS          WIDTH 80,60
' 1 TO 78 COLUMNS       WIDTH 80,60

PIXPERROW = 8: PIXPERCOL = 8

VUX% = 80: VUY% = 40: VLX% = 584: VLY% = 408       'VLX%=580
VIEW (VUX%, VUY%)-(VLX%, VLY%)
VXRANGE% = VLX% - VUX%: VYRANGE = VLY% - VUY%
LXRANGE = ABS(XMAX - XMIN): LYRANGE = ABS(YMAX - YMIN)
VYRANGE = VYRANGE * YF
XTOV = LXRANGE / VXRANGE%: YTOV = LYRANGE / VYRANGE
WXPIX% = 504: WYPIX% = 367
WX1 = X1 - 1 * 8 * XTOV: WY1 = Y1 + 16 * YTOV
WX2 = X2 + 1 * 8 * XTOV: WY2 = Y2 - 16 * YTOV
COLUMNSTART% = 11: COLUMNEND = 73: TOTALPHYSCOLS = 63

'  YROWSTART% = 7: YROWEND = 51  '45 ROWS
YROWSTART% = 7: YROWEND% = 56    '50 ROWS
WINDOW (WX1, WY1)-(WX2, WY2)

'LINE (WX1, WY1)-(WX2, WY2), GREEN, B   'GREEN BOX DEFINES PHYSICAL WINDOW
LINE (X1, Y1)-(X2, Y2), LIGHTBLUE%, B         'BLUE% BOX DEFINES LOGICAL WINDOW

COLOR RED%
LOCATE 1, 1: PRINT COMMONTITLE$
CALL SCENTER(84, NEWCOL, TITLE$): LOCATE 3, NEWCOL:  PRINT TITLE$
CALL SCENTER(84, NEWCOL, XAXISTITLE$): LOCATE 58, NEWCOL: PRINT XAXIS$
LOCATE 19, 1: PRINT YAXIS$
COLOR FORECOLOR

COLOR FORECOLOR%
XROW = 54: COUNT = 0
FOR I = X1 TO X2 STEP XSTEP  ' X AXIS TICK MARKS
LINE (I, Y2)-(I, WY2)
XCOL = 5 + I * LXPIXELTOCOL
IF I = X1 THEN : XCOL = 4 + I * LXPIXELTOCOL
'CALL NCENTER(XROW, XCOL, I)
CALL NPRINTFORMAT(I, TMP$)
CALL IFEVEN(COUNT, EVENORODD)  ' 2 IF EVEN, 1 IF ODD
IF EVENORODD = 2 THEN
XROW = 54
ELSE
XROW = 56
END IF
LOCATE XROW, XCOL: PRINT USING TMP$; I
COUNT = COUNT + 1
NEXT I

IF COUNT < 7 THEN 'BE SURE THAT XMAX PLOTS
I = X2: XROW = 54
LINE (I, Y2)-(I, WY2)
XCOL = 5 + I * LXPIXELTOCOL
CALL NPRINTFORMAT(I, TMP$)
LOCATE XROW, XCOL: PRINT USING TMP$; I
END IF

'COLOR GREEN%
COUNT = 0                       'Y AXIS
'IF YMIN >= 0 THEN               ' IF THERE IS NO ZERO FOR THE Y DATA
YDATAP = YMAX
FOR ROW = YROWSTART% TO YROWEND% STEP 5       ' 2
CALL NPRINTFORMAT(YDATAP, TMP$)
LOCATE ROW, 1:  PRINT USING TMP$; YDATAP
LINE (X1, YDATAP * YF)-(WX1, YDATAP * YF)
YDATAP = YDATAP - YSTEP  ' - (5 * YSTEP)
'IF YDATAP < Y1 THEN : EXIT FOR
COUNT = COUNT + 1
NEXT ROW
IF YDATAP < YMIN THEN
LINE (X1, YMIN * YF)-(WX1, YMIN * YF)
END IF
'END IF

IF COUNT < 9 THEN
YDATAP = Y2: ROW = 51
CALL NPRINTFORMAT(YDATAP, TMP$)
LOCATE ROW, 1:  PRINT USING TMP$; YDATAP
LINE (X1, YDATAP)-(WX1, YDATAP)
END IF

COLOR RED%     'ZERO LINES
LINE (WX1, 0)-(XMAX, 0), LIGHTBLUE%  'Y ZERO LINE
'LINE (0, YMIN * YF)-(0, YMAX * YF), GREEN   'X ZERO LINE

COLOR LIGHTGREEN%                        'PLOT THE DATA POINTS

XSTOP = DATATOTAL - 1
'XSTOP = XMAX
'FOR I = XMIN TO DATATOTAL
FOR I = LBOUND(X) TO XSTOP
'FOR I = XMIN TO XSTOP

IF X(I) > XMAX THEN : EXIT FOR
PSET (X(I) * XF, YF * Y(I))
IF PLOTPOINT% > 0 THEN
CIRCLE (X(I) * XF, YF * Y(I)), CIRCLERADIUS, MAGENTA
PAINT (X(I) * XF, YF * Y(I)), MAGENTA
END IF
NEXT I

IF PLOTPOINT% < 2 THEN
COLOR LIGHTRED%                         'CONNECT THE POINTS
FOR I = LBOUND(X) TO XSTOP - 1
IF X(I) > XMAX THEN : EXIT FOR
LINE (X(I), YF * Y(I))-(X(I + 1), YF * Y(I + 1))
NEXT I

ELSE  ' SPIKE THE DATA

FOR I = LBOUND(X) TO XMAX
IF X(I) > XMAX THEN : EXIT FOR
LINE (X(I), YF * Y(I))-(X(I), 0)
NEXT I
END IF

DUM$ = INPUT$(1)
VIEW: CLS 0
SCREEN 0
END SUB

REM $DYNAMIC
SUB GRAPHTEST

DO
CLS
DATACOUNT = 0: XDATATOTAL = 0: YDATATOTAL = 0

CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
CALL GETARRAYSIZE(INPUTFILE$, INFILENUM, XDATATOTAL, YDATATOTAL)
DATACOUNT = YDATATOTAL - 1
REDIM X(0 TO DATACOUNT), Y(0 TO DATACOUNT)


'INFILENUM = FREEFILE
OPEN INPUTFILE$ FOR INPUT AS #INFILENUM
IF XDATATOTAL = 1 THEN
FOR I = 0 TO DATACOUNT
INPUT #INFILENUM, Y(I)
X(I) = I
NEXT I
ELSE
FOR I = 0 TO DATACOUNT
INPUT #INFILENUM, X(I), Y(I)
NEXT I
END IF
CLOSE #INFILENUM

CALL GETMINMAX(X(), Y(), DATACOUNT, XMIN, XMAX, YMIN, YMAX)

PRINT "LBX "; LBOUND(X); " = "; X(LBOUND(X))
PRINT "UBX "; UBOUND(X); " = "; X(UBOUND(X))
PRINT "LBY "; LBOUND(Y); " = "; Y(LBOUND(Y))
PRINT "UBY "; UBOUND(Y); " = "; Y(UBOUND(Y))
'WAITFORKEY
CALL GRAPH(X(), Y(), DATACOUNT, XMIN, XMAX, YMIN, YMAX, 0, T$, X$, Y$)


INPUT "AGAIN "; YORN$
IF UCASE$(YORN$) <> "Y" THEN : EXIT DO
LOOP

END SUB

REM $STATIC
SUB GWHILBERT (HILBERTTRANSFORM(), W())
FOR I = LBOUND(HILBERTTRANSFORM) TO UBOUND(HILBERTTRANSFORM)
IF W(I) > 0 THEN
HILBERTTRANSFORM(I) = 2 * HILBERTTRANSFORM(I)
ELSEIF W(I) = 0 THEN
HILBERTTRANSFORM(I) = HILBERTTRANSFORM(I)
ELSE
HILBERTTRANSFORM(I) = 0  ' FOR W(I) < 0
END IF
NEXT I
END SUB

SUB HARM (A(), M(), INV(), S(), ifset, IFERR)
 'SUB HARM (a(), M%(), INV%(), S(), ifset%, IFERR%)
PRINT "STARTING HARM AT "; TIME$: INPUT DUM$
PRINT "IFSET IS "; ifset: INPUT DUM$

'DEFINT I-N
DIM NP(3), W(3), W2(3), W3(3)    'I ADDED THIS
10 IF (ABS(ifset) - 1) <= 1 THEN GOTO 900 ELSE GOTO 12
12        mtt = 0
        FOR J = 1 TO 3
        PRINT "M "; J; " IS "; M(J)
        IF M(J) > mtt THEN mtt = M(J)
        NEXT J
        mtt = mtt - 2
      ROOT2 = SQR(2!)
      IF mtt < 2 THEN GOTO 13 ELSE GOTO 14
      IF mtt > 12 THEN GOTO 13 ELSE GOTO 14
13 IFERR = 1
      PRINT "mtt is "; mtt
      PRINT "np too large or too small": INPUT DUM$
      GOTO 1000
14 IFERR = 0
      M1 = M(1)
      M2 = M(2)
      M3 = M(3)
      N1 = 2 ^ M1
      N2 = 2 ^ M2
      N3 = 2 ^ M3
16 IF ifset <= 0 THEN GOTO 18 ELSE GOTO 20
18 NX = N1 * N2 * N3
      QQN = NX
      FOR I = 1 TO NX
      A(2 * I - 1) = A(2 * I - 1) / QQN
19 A(2 * I) = -A(2 * I) / QQN
      NEXT I
20 NP(1) = N1 * 2
      NP(2) = NP(1) * N2
      NP(3) = NP(2) * N3
      FOR id = 1 TO 3       'end 250
      IL = NP(3) - NP(id)
      IL1 = IL + 1
      MI = M(id)
      IF MI <= 0 THEN GOTO 250
30 IDIF = NP(id)
      KBIT = NP(id)
      MEV = 2 * (MI / 2)
      IF (MI - MEV) <= 0 THEN GOTO 60
40 KBIT = KBIT / 2
      KL = KBIT - 2
      FOR I = 1 TO IL1 STEP IDIF  'end 50
      KLAST = KL + I
      FOR K = I TO KLAST STEP 2      'end 50
      KD = K + KBIT
      T = A(KD)
      A(KD) = A(K) - T
      A(K) = A(K) + T
      T = A(KD + 1)
      A(KD + 1) = A(K + 1) - T
50 A(K + 1) = A(K + 1) + T
      NEXT K
      NEXT I
      IF (MI - 1) <= 0 THEN GOTO 250 ELSE GOTO 52
52 LFIRST = 3
      JLAST = 1
      GOTO 70
60 LFIRST = 2
      JLAST = 0
70 FOR L = LFIRST TO MI STEP 2          'end 240
      JJDIF = KBIT
      KBIT = KBIT / 4
      KL = KBIT - 2
      FOR I = 1 TO IL1 STEP IDIF        'end 80
      KLAST = I + KL
      FOR K = I TO KLAST STEP 2            'end 80
      K1 = K + KBIT
      K2 = K1 + KBIT
      K3 = K2 + KBIT
      T = A(K2)
      A(K2) = A(K) - T
      A(K) = A(K) + T
      T = A(K2 + 1)
      A(K2 + 1) = A(K + 1) - T
      A(K + 1) = A(K + 1) + T
      T = A(K3)
      A(K3) = A(K1) - T
      A(K1) = A(K1) + T
      T = A(K3 + 1)
      A(K3 + 1) = A(K1 + 1) - T
      A(K1 + 1) = A(K1 + 1) + T
      T = A(K1)
      A(K1) = A(K) - T
      A(K) = A(K) + T
      T = A(K1 + 1)
      A(K1 + 1) = A(K + 1) - T
      A(K + 1) = A(K + 1) + T
      R = -A(K3 + 1)
      T = A(K3)
      A(K3) = A(K2) - R
      A(K2) = A(K2) + R
      A(K3 + 1) = A(K2 + 1) - T
80 A(K2 + 1) = A(K2 + 1) + T
      NEXT K
      NEXT I
      IF JLAST <= 0 THEN GOTO 235
82 JJ = JJDIF + 1
      ILAST = IL + JJ
      FOR I = JJ TO ILAST STEP IDIF            'end 85
      KLAST = KL + I
      FOR K = I TO KLAST STEP 2                 'end 85
      K1 = K + KBIT
      K2 = K1 + KBIT
      K3 = K2 + KBIT
      R = -A(K2 + 1)
      T = A(K2)
      A(K2) = A(K) - R
      A(K) = A(K) + R
      A(K2 + 1) = A(K + 1) - T
      A(K + 1) = A(K + 1) + T
      AWR = A(K1) - A(K1 + 1)
      AWI = A(K1 + 1) + A(K1)
      R = -A(K3) - A(K3 + 1)
      T = A(K3) - A(K3 + 1)
      A(K3) = (AWR - R) / ROOT2
      A(K3 + 1) = (AWI - T) / ROOT2
      A(K1) = (AWR + R) / ROOT2
      A(K1 + 1) = (AWI + T) / ROOT2
      T = A(K1)
      A(K1) = A(K) - T
      A(K) = A(K) + T
      T = A(K1 + 1)
      A(K1 + 1) = A(K + 1) - T
      A(K + 1) = A(K + 1) + T
      R = -A(K3 + 1)
      T = A(K3)
      A(K3) = A(K2) - R
      A(K2) = A(K2) + R
      A(K3 + 1) = A(K2 + 1) - T
85 A(K2 + 1) = A(K2 + 1) + T
      NEXT K
      NEXT I
      IF (JLAST - 1) <= 0 THEN GOTO 235 ELSE GOTO 90
90 JJ = JJ + JJDIF
      FOR J = 2 TO JLAST           'end 230
96 I = INV(J + 1)
98 IC = NT - I
      W(1) = S(IC)
      W(2) = S(I)
      I2 = 2 * I
      i2c = NT - I2
      IF i2c < 0 THEN GOTO 120
      IF i2c = 0 THEN GOTO 110 ELSE GOTO 100
100 W2(1) = S(i2c)
      W2(2) = S(I2)
      GOTO 130
110 W2(1) = 0!
      W2(2) = 1!
      GOTO 130
120 I2CC = i2c + NT
      i2c = -i2c
      W2(1) = -S(i2c)
      W2(2) = S(I2CC)
130 I3 = I + I2
      I3C = NT - I3
      IF I3C < 0 THEN GOTO 160
      IF I3C = 0 THEN GOTO 150 ELSE GOTO 140
140 W3(1) = S(I3C)
      W3(2) = S(I3)
      GOTO 200
150 W3(1) = 0!
      W3(2) = 1!
      GOTO 200
160 i3cc = I3C + NT
      IF i3cc < 0 THEN GOTO 190
      IF i3cc = 0 THEN GOTO 180 ELSE 170
170 I3C = -I3C
      W3(1) = -S(I3C)
      W3(2) = S(i3cc)
      GOTO 200
180 W3(1) = -1!
      W3(2) = 0!
      GOTO 200
190 I3CCC = NT + i3cc
      i3cc = -i3cc
      W3(1) = -S(I3CCC)
      W3(2) = -S(i3cc)
200 ILAST = IL + JJ
      FOR I = JJ TO ILAST STEP IDIF       'end 220
      KLAST = KL + I
      FOR K = I TO KLAST STEP 2          'end 220
      K1 = K + KBIT
      K2 = K1 + KBIT
      K3 = K2 + KBIT
      R = A(K2) * W2(1) - A(K2 + 1) * W2(2)
      T = A(K2) * W2(2) + A(K2 + 1) * W2(1)
      A(K2) = A(K) - R
      A(K) = A(K) + R
      A(K2 + 1) = A(K + 1) - T
      A(K + 1) = A(K + 1) + T
      R = A(K3) * W3(1) - A(K3 + 1) * W3(2)
      T = A(K3) * W3(2) + A(K3 + 1) * W3(1)
      AWR = A(K1) * W(1) - A(K1 + 1) * W(2)
      AWI = A(K1) * W(2) + A(K1 + 1) * W(1)
      A(K3) = AWR - R
      A(K3 + 1) = AWI - T
      A(K1) = AWR + R
      A(K1 + 1) = AWI + T
      T = A(K1)
      A(K1) = A(K) - T
      A(K) = A(K) + T
      T = A(K1 + 1)
      A(K1 + 1) = A(K + 1) - T
      A(K + 1) = A(K + 1) + T
      R = -A(K3 + 1)
      T = A(K3)
      A(K3) = A(K2) - R
      A(K2) = A(K2) + R
      A(K3 + 1) = A(K2 + 1) - T
220 A(K2 + 1) = A(K2 + 1) + T
       NEXT K
       NEXT I
230 JJ = JJDIF + JJ
        NEXT J
235 JLAST = 4 * JLAST + 3
REM   240 CONTINUE
      NEXT L
250   NEXT id
REM   250 CONTINUE
      NTSQ = NT * NT
      M3MT = M3 - mt
350 IF M3MT < 0 THEN GOTO 370
360 igo3 = 1
      N3VNT = N3 / NT
      MINN3 = NT
      GOTO 380
370 igo3 = 2
      N3VNT = 1
      NTVN3 = NT / N3
      MINN3 = N3
380 JJD3 = NTSQ / N3
      M2MT = M2 - mt
450 IF M2MT < 0 THEN GOTO 470
460 igo2 = 1
      N2VNT = N2 / NT
      MINN2 = NT
      GOTO 480
470 igo2 = 2
      N2VNT = 1
      NTVN2 = NT / N2
      MINN2 = N2
480 JJD2 = NTSQ / N2
      M1MT = M1 - mt
550 IF M1MT < 0 THEN GOTO 570
560 igo1 = 1
      N1VNT = N1 / NT
      MINN1 = NT
      GOTO 580
570 igo1 = 2
      N1VNT = 1
      NTVN1 = NT / N1
      MINN1 = N1
580 JJD1 = NTSQ / N1
600 JJ3 = 1
      J = 1
      FOR jpp3 = 1 TO N3VNT     'end 880
      IPP3 = INV(JJ3)
      FOR jp3 = 1 TO MINN3     'end 870
      IF igo3 = 2 THEN GOTO 620        'GO TO (610,620),IGO3
610 IP3 = INV(jp3) * N3VNT
      GOTO 630
620 IP3 = INV(jp3) / NTVN3
630 I3 = (IPP3 + IP3) * N2
700 JJ2 = 1
      FOR jpp2 = 1 TO N2VNT       'end 870
      IPP2 = INV(JJ2) + I3
      FOR jp2 = 1 TO MINN2        'end 860
      IF igo2 = 2 THEN GOTO 720    'GO TO (710,720),IGO2
710 IP2 = INV(jp2) * N2VNT
      GOTO 730
720 IP2 = INV(jp2) / NTVN2
730 I2 = (IPP2 + IP2) * N1
800 JJ1 = 1
      FOR jpp1 = 1 TO N1VNT           'end 860
      IPP1 = INV(JJ1) + I2
      FOR jp1 = 1 TO MINN1            'end 850
      IF igo1 = 2 THEN GOTO 820 ELSE GOTO 810  'GO TO (810,820),IGO1
810 IP1 = INV(jp1) * N1VNT
      GOTO 830
820 IP1 = INV(jp1) / NTVN1
830 I = 2 * (IPP1 + IP1) + 1
      IF (J - I) < 0 THEN GOTO 840 ELSE GOTO 845
840 T = A(I)
      A(I) = A(J)
      A(J) = T
      T = A(I + 1)
      A(I + 1) = A(J + 1)
      A(J + 1) = T
845 AQAQ = 2
850 J = J + 2
    NEXT jp1
860 JJ1 = JJ1 + JJD1
    NEXT jpp1
    NEXT jp2
870 JJ2 = JJ2 + JJD2
    NEXT jpp2
    NEXT jp3
880 JJ3 = JJ3 + JJD3
    NEXT jpp3
890 IF ifset < 0 THEN GOTO 891 ELSE GOTO 895
891 FOR I = 1 TO NX               'end 892
892 A(2 * I) = -A(2 * I)
    NEXT I
895 GOTO 1000
900 mt = M(1) - 2: IFERR = 0
      NT = 2 ^ mt
      NTV2 = NT / 2
910 THETA = .7853981634#
      JSTEP = NT
      JDIF = NTV2
      S(JDIF) = SIN(THETA)
      FOR L = 2 TO mt     'end 950
      THETA = THETA / 2!
      JSTEP2 = JSTEP
      JSTEP = JDIF
      JDIF = JSTEP / 2
      S(JDIF) = SIN(THETA)
      JC1 = NT - JDIF
      S(JC1) = COS(THETA)
      JLAST = NT - JSTEP2
      IF (JLAST - JSTEP) < 0 THEN GOTO 950
920 FOR J = JSTEP TO JLAST STEP JSTEP        'end 940
      JC = NT - J
      JD = J + JDIF
940 S(JD) = S(J) * S(JC1) + S(JDIF) * S(JC)
     NEXT J
REM  950 CONTINUE
950  NEXT L
960  MTLEXP = NTV2
      LM1EXP = 1
      INV(1) = 0
      FOR L = 1 TO mt    'end 980
      INV(LM1EXP + 1) = MTLEXP
      FOR J = 2 TO LM1EXP      'end 970
      JJ = J + LM1EXP
970 INV(JJ) = INV(J) + MTLEXP
      NEXT J
      MTLEXP = MTLEXP / 2
980 LM1EXP = LM1EXP * 2
     NEXT L
982 IF ifset = 0 THEN GOTO 895 ELSE GOTO 12
1000 AQAQ = 3

PRINT "ENDING HARM AT "; TIME$: INPUT DUM$
END SUB

SUB IFEVEN (NUM1, EVENORODD) ' 2 IF EVEN, 1 IF ODD
NUM2 = NUM1 MOD 2
IF NUM2 = 0 THEN : EVENORODD = 2: EXIT SUB
EVENORODD = 1
END SUB

SUB LEASTDELTA (XGET(), LEAST)
LEAST = XGET(UBOUND(XGET)) - XGET(LBOUND(XGET))
FOR I = (LBOUND(XGET)) TO (UBOUND(XGET) - 1)
TEST = (XGET(I + 1) - XGET(I))
IF TEST < LEAST THEN LEAST = TEST
NEXT I
END SUB

FUNCTION LN (X)
IF X > 0 THEN
LN = LOG(X)
ELSE
PRINT "CANNOT TAKE LN(X<0)"
EXIT FUNCTION
END IF
END FUNCTION

SUB MATCHFILTER
'NOTE:  SINCE THE CONVOLUTION ROUTINE I'M USING FOLDS THE 2nd SIGNAL, I
'       NEED TO FOLD THE 2nd SIGNAL IN THE CORRELATION ROUTINE BEFORE
'       I SEND IT TO THE CONVOLUTION IN ORDER TO REVERSE THE EFFECTS OF
'       THE FOLDING IN THE CONVOLUTION ! BUT, WITH THE MATCH FILTER, WE
'       WANT TO FOLD THE 2nd SIGNAL!  SO IN THIS MATCH FILTER ROUTINE,
'       THE FOLDING IS DONE ONLY FOR THE PURPOSE OF DISPLAYING THE
'       GRAPH OF THE FOLD !

SNP = 0: GNP = 0: ISNP = 0: IGNP = 0: CONTOTAL = 0: DIVISOR = 1
PRINT "MATCH FILTER"
PRINT "SOLVE FOR h(t) = s(t) * g(t)"
PRINT "FOR THE SIGNAL S(t)"
CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
SINPUTFILE$ = INPUTFILE$: SINFILENUM = INFILENUM
CALL GETOUTPUTFILE(OUTPUTFILE$, OUTFILENUM)
'INPUTFILE$ = "S.DAT"
CALL GETARRAYSIZE(SINPUTFILE$, SINFILENUM, XDATATOTAL, YDATATOTAL)
ISNP = YDATATOTAL - 1
PRINT "ISNP,YT "; ISNP, YDATATOTAL

PRINT "FOR THE SIGNAL G(t-T)"
CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
GINPUTFILE$ = INPUTFILE$: GINFILENUM = INFILENUM
CALL GETARRAYSIZE(GINPUTFILE$, GINFILENUM, XDATATOTAL, YDATATOTAL)
IGNP = YDATATOTAL - 1
PRINT "IGNP, YT "; IGNP, YDATATOTAL

'**************************************

IF ISNP < IGNP THEN
GNP = IGNP
SNP = GNP
ELSEIF IGNP < ISNP THEN
SNP = ISNP
GNP = SNP
ELSE
SNP = ISNP
GNP = SNP
END IF
PRINT "ISNP "; ISNP; " IGNP "; IGNP; " SNP "; SNP; " GNP "; GNP

START% = 0
REDIM S(START% TO SNP), SCOUNT(START% TO SNP), CON(START% TO SNP * 2)
REDIM G(START% TO GNP), GCOUNT(START% TO GNP), GFOLD(START% TO GNP)

FOR I = 0 TO SNP - 1
SCOUNT(I) = I
GCOUNT(I) = I
NEXT I

IF ISNP < IGNP THEN
FOR I = ISNP + 1 TO GNP
S(I) = 0
NEXT I
ELSEIF IGNP < ISNP THEN
FOR I = IGNP + 1 TO SNP
G(I) = 0
NEXT I
END IF

PRINT "ISNP "; ISNP; " SNP "; SNP
PRINT "IGNP "; IGNP; " GNP "; GNP

OPEN SINPUTFILE$ FOR INPUT AS #SINFILENUM
FOR I = 0 TO ISNP - 1
INPUT #SINFILENUM, S(I)    'THE SIGNAL
NEXT I
CLOSE #SINFILENUM

OPEN GINPUTFILE$ FOR INPUT AS #GINFILENUM
FOR I = 0 TO IGNP - 1
INPUT #GINFILENUM, G(I)    'THE SIGNAL
NEXT I
CLOSE #GINFILENUM

COLOR RED%: PRINT "WORKING...": COLOR GREEN%

CALL FOLDANDSHIFT(GNP, G(), GFOLD())

'PRINT "I", "S", "G", "GF"
'FOR I = 0 TO SNP - 1
'PRINT I, S(I), G(I), GFOLD(I)
'NEXT I

'*****************************************
LX = SNP: LB = GNP
LY = 0: LY = LX + LB: CONTOTAL = LY
REDIM LX(1 TO LX), LB(1 TO LB), CY(1 TO LY)

PRINT "LX "; LX; " LB "; LB; " LY "; LY
REDIM CONCOUNT(0 TO CONTOTAL)
FOR I = 0 TO CONTOTAL - 1
CONCOUNT(I) = I
NEXT I
FOR I = 1 TO SNP
LX(I) = S(I - 1)
LB(I) = G(I - 1)
NEXT I

FOR I = 1 TO LX
FOR J = 1 TO LB
CY(I + J - 1) = CY(I + J - 1) + LX(I) * LB(J)
NEXT J
NEXT I
PRINT "IN NEWCONV"
'FOR I = 1 TO 2 * SNP
'PRINT CY(I)
'NEXT I
FOR I = 0 TO SNP * 2 - 1
CON(I) = CY(I + 1)
NEXT I
ERASE CY, LX, LB
'*****************************************

PRINT "ISNP "; ISNP; " SNP "; SNP
PRINT "IGNP "; IGNP; " GNP "; GNP
PRINT "CONTOTAL "; CONTOTAL
INPUT DUM$

PRINT "OUTPUT WRITTEN TO "; OUTPUTFILE$
OPEN OUTPUTFILE$ FOR OUTPUT AS #OUTFILENUM

PRINT #OUTFILENUM, "I", , "S(I)", "G(I)"
TMP$ = "############.#######"
FOR I = 1 TO ISNP
PRINT #OUTFILENUM, I,
PRINT #OUTFILENUM, USING TMP$; S(I); G(I)
NEXT I
PRINT #OUTFILENUM, "I", "CON(I)"
FOR I = 1 TO CONTOTAL
PRINT #OUTFILENUM, I,
PRINT #OUTFILENUM, USING TMP$; CON(I)
NEXT I
CLOSE

PRINT "READY TO GRAPH S(t)"
P = SNP
CALL GETMINMAX(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX)
'INPUT DUM$
TITLE$ = "S(t)"
YTITLE$ = "S(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH g(t)"
P = GNP
CALL GETMINMAX(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "g(t)"
YTITLE$ = "g(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE FOLDED SIGNAL"
P = GNP
CALL GETMINMAX(GCOUNT(), GFOLD(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "G(t) FOLDED"
YTITLE$ = "G(t) FOLDED"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), GFOLD(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)


PRINT "READY TO GRAPH THE CONVOLUTION MATCH FILTER"
'LINE INPUT "ENTER THE TITLE "; TITLE$
TITLE$ = "MATCH FILTER OF "
INPUT "MATCH FILTER OF "; TCC$
TITLE$ = TITLE$ + TCC$

CALL GETMINMAX(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX)
PRINT "CONTOTAL "; CONTOTAL
PRINT "XMIN "; XMIN
PRINT "XMAX "; XMAX
'TITLE$ = "CONVOLUTION p(t)=s(t) * g(t)"
YTITLE$ = "P(t)"
XTITLE$ = ""
CALL GRAPH(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

'WAITFORKEY
ERASE S, SCOUNT, CON, G, GCOUNT

END SUB

SUB MENUSIGNALPRO
CLS : DO: COLOR GREEN%: PRINT

PRINT : PRINT " 1) AMPLITUDE SPECTRUM"
PRINT : PRINT " 2) FAST FOURIER TRANSFORM"
PRINT : PRINT " 3) BAND PASS FILTER"
PRINT : PRINT " 4) CONVOLUTION"
PRINT : PRINT " 5) LOW PASS FILTER BY CONVOLUTION"
PRINT : PRINT " 6) CROSS CORRELATION"
PRINT : PRINT " 7) MATCH FILTER"
PRINT : PRINT " 8) DIMENSIONAL TRANSFORMS"
PRINT : PRINT " 9) DE-GHOSTING"
PRINT : PRINT "10) QUIT THIS MENU"
PRINT : INPUT "ENTER A NUMBER TO INDICATE YOUR CHOICE "; CHOICE%

SELECT CASE CHOICE%
        CASE 1:


        CASE 2:
        CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
        CALL GETOUTPUTFILE(OUTPUTFILE$, OUTFILENUM)
        CALL GETARRAYSIZE(INPUTFILE$, INFILENUM, XDATATOTAL, YDATATOTAL)
        NP = YDATATOTAL - 1
        CALL FASTFOURIER(INPUTFILE$, INFILENUM, OUTPUTFILE$, OUTFILENUM, NP)

        CASE 3:

'CALL BANDPASSFILTER(S(), TIME(), CN(), FREQ(), RA(), IA(), PHASEANGLE(), NP2)

        CASE 4:   CALL CONVOLUTION

        CASE 5:    'FILTER BY CONVOLUTION

        CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
        CALL GETOUTPUTFILE(OUTPUTFILE$, OUTFILENUM)
        CALL GETARRAYSIZE(INPUTFILE$, INFILENUM, XDATATOTAL, YDATATOTAL)
        NP = YDATATOTAL - 1
INPUT "ENTER THE SAMPLE RATE IN MSEC (1.0) "; R$: CALL NIFDEFAULT(R$, DELTAT, 1)
        DELTAT = DELTAT * .001


                        ISNP = NP
                        REDIM S(0 TO ISNP - 1), FILT(0 TO ISNP - 1)
                        REDIM FOLD(0 TO ISNP - 1), CON(0 TO ISNP * 2)
                        REDIM CONCOUNT(0 TO ISNP * 2)
                        INFILENUM = FREEFILE
                        OPEN INPUTFILE$ FOR INPUT AS #INFILENUM
                        FOR I = 0 TO ISNP - 1
                        INPUT #INFILENUM, S(I)
                        NEXT I
                        CLOSE #INFILENUM

                        CALL SINCLPF(ISNP, DELTAT, FILT())

                        FOR I = ISNP + 1 TO NP2 - 1
                        FILT(I) = 0
                        S(I) = 0
                        NEXT I

                CALL FOLDANDSHIFT(ISNP, FILT(), FOLD())

                CALL NEWCONV(ISNP, ISNP, S(), FILT(), FOLD(), CON(), CONCOUNT())

        CASE 6:
                CALL CROSSCORRELATION



        CASE 7:
                CALL MATCHFILTER

        CASE 8:
                 CALL DIMTRANSFORM

        CASE 9:
               ' COLOR 9
               ' PRINT "TYPE EXIT TO RETURN TO PROGRAM"
               ' SHELL
               ' COLOR 2
                CALL DEGHOST


        CASE 10: DONE = TRUE


END SELECT       'SELECT CASE CHOICE
LOOP UNTIL DONE


END SUB

SUB NCENTER (ROW, NEWCOL, NUMBERLABLE)
'x% = x% + Down%
NEWCOL = INT((NEWCOL - LEN(NUMBERLABLE)) / 2)
'LOCATE ROW, COLUMN: PRINT NUMBERLABLE
END SUB

REM $DYNAMIC
SUB NEWCONV (ISNP, IGNP, S(), G(), GF(), CON(), CONCOUNT())
'THE ARRAYS S(), G(), AND GF() MUST BE THE SAME SIZE BEFORE THEY GET HERE
IF ISNP < IGNP THEN
GNP = IGNP
SNP = GNP
PRINT "ISNP < IGNP"
ELSEIF IGNP < ISNP THEN
SNP = ISNP
GNP = SNP
PRINT "IGNP < ISNP"
ELSE
PRINT "ISNP=IGNP"
SNP = ISNP
GNP = SNP
END IF

PRINT "ISNP "; ISNP; " IGNP "; IGNP; " SNP "; SNP; " GNP "; GNP

REDIM GCOUNT(0 TO GNP - 1), SCOUNT(0 TO SNP - 1)

IF ISNP < IGNP THEN
FOR I = ISNP + 1 TO GNP - 1
S(I) = 0
NEXT I
ELSEIF IGNP < ISNP THEN
FOR I = IGNP + 1 TO SNP - 1
G(I) = 0
NEXT I
END IF

PRINT "ISNP "; ISNP; " SNP "; SNP
PRINT "IGNP "; IGNP; " GNP "; GNP

FOR I = 0 TO SNP - 1
SCOUNT(I) = I
GCOUNT(I) = I
NEXT I

'*****************************************
LX = SNP: LB = GNP
LY = 0: LY = LX + LB: CONTOTAL = LY
REDIM LX(1 TO LX), LB(1 TO LB), CY(1 TO LY)

PRINT "LX "; LX; " LB "; LB; " LY "; LY
PRINT "CONTOTAL "; CONTOTAL
'INPUT DUM$
REDIM CONCOUNT(0 TO CONTOTAL)
FOR I = 0 TO CONTOTAL
CONCOUNT(I) = I
NEXT I
FOR I = 1 TO SNP
LX(I) = S(I - 1)
LB(I) = G(I - 1)
NEXT I

TMP$ = "###.###": STARTTIME$ = TIME$
COUNT = 0
FOR I = 1 TO LX
FOR J = 1 TO LB
CY(I + J - 1) = CY(I + J - 1) + LX(I) * LB(J)
COUNT = COUNT + 1
NEXT J
CLS
PRINT "STARTED CONVOLUTION AT "; STARTTIME$
PRINT USING TMP$; COUNT / (SNP * .5 * CONTOTAL) * 100;
PRINT " % COMPLETED"
NEXT I
PRINT "ENDED CONVOLUTION AT "; TIME$
PRINT "NUMBER OF CONVOLUTION CALCULATIONS "; COUNT
FOR I = 0 TO SNP * 2 - 1
CON(I) = CY(I + 1)
NEXT I
ERASE CY, LX, LB

'PRINT "CONTOTAL "; CONTOTAL: INPUT DUM$
'*****************************************
'PRINT "READY TO GRAPH S(t)"
'P = SNP
'CALL GETMINMAX(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX)
'TITLE$ = "S(t)"
'YTITLE$ = "S(t)"
'XTITLE$ = "POINT COUNT"
'CALL GRAPH(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$,TSCOUNT,)

'PRINT "READY TO GRAPH g(t)"
'P = GNP
'CALL GETMINMAX(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX)
'TITLE$ = "g(t)"
'YTITLE$ = "g(t)"
'XTITLE$ = "POINT COUNT"
'CALL GRAPH(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$,TSCOUNT,)

'PRINT "READY TO GRAPH THE FOLDED SIGNAL"
'P = GNP
'CALL GETMINMAX(GCOUNT(), GF(), P, XMIN, XMAX, YMIN, YMAX)
'TITLE$ = "G(t) FOLDED"
'YTITLE$ = "G(t) FOLDED"
'XTITLE$ = "POINT COUNT"
'CALL GRAPH(GCOUNT(), GF(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$,TSCOUNT,)

'PRINT "READY TO GRAPH THE CONVOLUTION"
'LINE INPUT "ENTER THE TITLE "; TITLE$
'CALL GETMINMAX(CONCOUNT(), CON(), CONTOTAL - 1, XMIN, XMAX, YMIN, YMAX)
'PRINT "CONTOTAL "; CONTOTAL
'PRINT "XMIN "; XMIN
'PRINT "XMAX "; XMAX
'TITLE$ = "CONVOLUTION p(t)=s(t) * g(t)"
'YTITLE$ = "P(t)"
'XTITLE$ = ""
'CALL GRAPH(CONCOUNT(), CON(), CONTOTAL - 1, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$,TSCOUNT,)

ERASE SCOUNT, GCOUNT
END SUB

REM $STATIC
SUB NEWCONVOLUTION
SNP = 0: GNP = 0: ISNP = 0: IGNP = 0: CONTOTAL = 0: DIVISOR = 1
PRINT "SOLVE FOR h(t) = s(t) * g(t)"
PRINT "FOR THE SIGNAL S(t)"
CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
SINPUTFILE$ = INPUTFILE$: SINFILENUM = INFILENUM
CALL GETOUTPUTFILE(OUTPUTFILE$, OUTFILENUM)
'INPUTFILE$ = "S.DAT"
CALL GETARRAYSIZE(SINPUTFILE$, SINFILENUM, XDATATOTAL, YDATATOTAL)
ISNP = YDATATOTAL - 1
PRINT "ISNP,YT "; ISNP, YDATATOTAL

PRINT "FOR THE SIGNAL G(t-T)"
CALL GETINPUTFILE(INPUTFILE$, INFILENUM)
GINPUTFILE$ = INPUTFILE$: GINFILENUM = INFILENUM
CALL GETARRAYSIZE(GINPUTFILE$, GINFILENUM, XDATATOTAL, YDATATOTAL)
IGNP = YDATATOTAL - 1
PRINT "IGNP, YT "; IGNP, YDATATOTAL

'**************************************

IF ISNP < IGNP THEN
GNP = IGNP
SNP = GNP
ELSEIF IGNP < ISNP THEN
SNP = ISNP
GNP = SNP
ELSE
SNP = ISNP
GNP = SNP
END IF
PRINT "ISNP "; ISNP; " IGNP "; IGNP; " SNP "; SNP; " GNP "; GNP

START% = 0
REDIM S(START% TO SNP), SCOUNT(START% TO SNP), CON(START% TO SNP * 2)
REDIM G(START% TO GNP), GCOUNT(START% TO GNP), GFOLD(START% TO GNP)

FOR I = 0 TO SNP - 1
SCOUNT(I) = I
GCOUNT(I) = I
NEXT I

IF ISNP < IGNP THEN
FOR I = ISNP + 1 TO GNP
S(I) = 0
NEXT I
ELSEIF IGNP < ISNP THEN
FOR I = IGNP + 1 TO SNP
G(I) = 0
NEXT I
END IF

PRINT "ISNP "; ISNP; " SNP "; SNP
PRINT "IGNP "; IGNP; " GNP "; GNP

OPEN SINPUTFILE$ FOR INPUT AS #SINFILENUM
FOR I = 0 TO ISNP - 1
INPUT #SINFILENUM, S(I)    'THE SIGNAL
NEXT I
CLOSE #SINFILENUM

OPEN GINPUTFILE$ FOR INPUT AS #GINFILENUM
FOR I = 0 TO IGNP - 1
INPUT #GINFILENUM, G(I)    'THE SIGNAL
NEXT I
CLOSE #GINFILENUM

COLOR RED%: PRINT "WORKING...": COLOR GREEN%
CALL FOLDANDSHIFT(GNP, G(), GFOLD())

PRINT "I", "S", "G", "GF"
FOR I = 0 TO SNP - 1
PRINT I, S(I), G(I), GFOLD(I)
NEXT I

'*****************************************
LX = SNP: LB = GNP
LY = 0: LY = LX + LB: CONTOTAL = LY
REDIM LX(1 TO LX), LB(1 TO LB), CY(1 TO LY)

PRINT "LX "; LX; " LB "; LB; " LY "; LY
REDIM CONCOUNT(0 TO CONTOTAL)
FOR I = 0 TO CONTOTAL - 1
CONCOUNT(I) = I
NEXT I
FOR I = 1 TO SNP
LX(I) = S(I - 1)
LB(I) = G(I - 1)
NEXT I

FOR I = 1 TO LX
FOR J = 1 TO LB
CY(I + J - 1) = CY(I + J - 1) + LX(I) * LB(J)
NEXT J
NEXT I
PRINT "IN NEWCONV"
FOR I = 1 TO 2 * SNP
PRINT CY(I)
NEXT I
FOR I = 0 TO SNP * 2 - 1
CON(I) = CY(I + 1)
NEXT I
ERASE CY, LX, LB
'*****************************************

PRINT "ISNP "; ISNP; " SNP "; SNP
PRINT "IGNP "; IGNP; " GNP "; GNP
PRINT "CONTOTAL "; CONTOTAL
INPUT DUM$

PRINT "OUTPUT WRITTEN TO "; OUTPUTFILE$
OPEN OUTPUTFILE$ FOR OUTPUT AS #OUTFILENUM

PRINT #OUTFILENUM, "I", , "S(I)", "G(I)"
TMP$ = "############.#######"
FOR I = 1 TO ISNP
PRINT #OUTFILENUM, I,
PRINT #OUTFILENUM, USING TMP$; S(I); G(I)
NEXT I
PRINT #OUTFILENUM, "I", "CON(I)"
FOR I = 1 TO CONTOTAL
PRINT #OUTFILENUM, I,
PRINT #OUTFILENUM, USING TMP$; CON(I)
NEXT I
CLOSE

PRINT "READY TO GRAPH S(t)"
P = SNP
CALL GETMINMAX(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX)
'INPUT DUM$
TITLE$ = "S(t)"
YTITLE$ = "S(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(SCOUNT(), S(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH g(t)"
P = GNP
CALL GETMINMAX(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "g(t)"
YTITLE$ = "g(t)"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), G(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE FOLDED SIGNAL"
P = GNP
CALL GETMINMAX(GCOUNT(), GFOLD(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "G(t) FOLDED"
YTITLE$ = "G(t) FOLDED"
XTITLE$ = "POINT COUNT"
CALL GRAPH(GCOUNT(), GFOLD(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)


PRINT "READY TO GRAPH THE CONVOLUTION"
'LINE INPUT "ENTER THE TITLE "; TITLE$
TITLE$ = "Trace1 * f(t), f0 = 40 Hz"
CALL GETMINMAX(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX)
PRINT "CONTOTAL "; CONTOTAL
PRINT "XMIN "; XMIN
PRINT "XMAX "; XMAX
'TITLE$ = "CONVOLUTION p(t)=s(t) * g(t)"
YTITLE$ = "P(t)"
XTITLE$ = ""
CALL GRAPH(CONCOUNT(), CON(), CONTOTAL, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)


'WAITFORKEY
ERASE S, SCOUNT, CON, G, GCOUNT
END SUB

SUB NEWFOLD (IP, SIG(), FOLD())
M = 1
FOR J = IP TO 1 STEP -1
FOLD(M) = SIG(J)
M = M + 1
NEXT J



END SUB

SUB NIFDEFAULT (R$, R, DEFAULTVAL)
'CALL NIFDEFAULT(R$, DELTAT, 1)

IF R$ = "" THEN
R = DEFAULTVAL: EXIT SUB
ELSE
R = VAL(R$)
END IF
END SUB

SUB NPRINTFORMAT (NUM, PRNFMT$)
NUM$ = STR$(NUM)
TEST = ABS(NUM)
IF TEST = 0 THEN
PRNFMT$ = "#########": EXIT SUB
'ELSEIF NUM < .00001 THEN
'NUM = 0
'PRNFMT$ = "#": EXIT SUB
ELSEIF TEST < .0001 THEN
PRNFMT$ = "##.##^^^^": EXIT SUB
ELSEIF TEST < 1 THEN
PRNFMT$ = "##.######": EXIT SUB
ELSEIF TEST < 10000 THEN
PRNFMT$ = "#####.###": EXIT SUB
ELSEIF TEST < 100000 THEN
PRNFMT$ = "######.##": EXIT SUB
ELSEIF TEST < 1000000 THEN
PRNFMT$ = "##.##^^^^": EXIT SUB
ELSE
PRNFMT$ = "##.#^^^^^": EXIT SUB
END IF

END SUB

SUB PAGE STATIC

PAGECOUNT = PAGECOUNT + 1
IF PAGECOUNT = 15 THEN
WAITFORKEY
PAGECOUNT = 0
CLS
END IF

END SUB

SUB PHASE (A, B, PHI)

IF B = 0 THEN
        IF A = 0 THEN PHI = 0: EXIT SUB
        IF A >= 0 THEN PHI = PI / 2: EXIT SUB
        IF A < 0 THEN PHI = (3 / 2) * PI: EXIT SUB
ELSEIF A >= 0 THEN
        IF B > 0 THEN PHI = ATN(A / B): EXIT SUB
        IF B < 0 THEN PHI = PI + ATN(A / B): EXIT SUB
ELSEIF A < 0 THEN
        IF B > 0 THEN PHI = 2 * PI + ATN(A / B): EXIT SUB
        IF B < 0 THEN PHI = PI + ATN(A / B): EXIT SUB
END IF

'IF (A = 0 AND B = 0) THEN PHI = 0: EXIT SUB
'IF A > 0 AND B = 0 THEN PHI = PI / 2: EXIT SUB
'IF A < 0 AND B = 0 THEN PHI = (3 / 2 * PI): EXIT SUB

'IF A >= 0 AND B > 0 THEN              'QUAD 1
'PHI = ATN(A / B): EXIT SUB

'ELSEIF A >= 0 AND B = 0 THEN
'PHI = PI / 2

'ELSEIF A < 0 AND B > 0 THEN          'QUAD 2
'PHI = 2 * PI + ATN(A / B): EXIT SUB

'ELSEIF A < 0 AND B < 0 THEN          'QUAD 3
'PHI = PI + ATN(A / B): EXIT SUB

'ELSEIF (A > 0) AND (B < 0) THEN           'QUAD 4
'PHI = PI + ATN(A / B): EXIT SUB

'END IF


END SUB

SUB PRINTCENTERED (S$, WIDE%)
L% = LEN(S$)
IF L% >= WIDE% THEN : EXIT SUB


S$ = SPACE$((WIDE% - L%) / 2) + S$


END SUB

SUB RESETALL

ROWLENGTH = 0
ROWLENGTH$ = "0"
DATACOUNT = 0
XDATATOTAL = 0
YDATATOTAL = 0
COMMAS = 0
ROWS = 0
COLUMNS = 0
ROWTOTAL = 0
COLUMNTOTAL = 0
XCOUNT = 0
YCOUNT = 0
XMIN = 0
XMAX = 0
YMIN = 0
YMAX = 0
SUM = 0
TMP = 0


END SUB

SUB SCENTER (OLDCOL, NEWCOL, LABLE$)
NEWCOL = INT((OLDCOL - LEN(LABLE$)) / 2)
END SUB

SUB SCREENTEST


SCREEN 12
WIDTH 80, 60
VIEW (1, 1)-(639, 479)

FOR ROW = 1 TO 58  '58 ROWS TOTAL
LOCATE ROW, 1: PRINT ROW
NEXT ROW

FOR COL = 3 TO 78 STEP 15     '78 COLS
LOCATE 58, COL: PRINT 0
NEXT COL

'LOCATE 2, 35: PRINT "TITLE"
'LOCATE 28, 35: PRINT "X AXIS"
'LOCATE 12, 1: PRINT "Y"




WAITFORKEY
SCREEN 0

END SUB

SUB SEISMO2 (NEWMODEL)
'***********************  GET E(t)
DEL$ = CHR$(127)
PIE$ = CHR$(227)
PHI$ = CHR$(237)
RHO$ = CHR$(235)

NLAYERS = 0:

PRINT "PHI "; PHI$; " DEL "; DEL$; " PI "; PIE$

IF NEWMODEL THEN
'INPUT "ENTER DELTA T IN SECONDS "; DELTAT
INPUT "ENTER THE NUMBER OF LAYERS INCLUDING THE HALF SPACE "; NLAYERS
REDIM V(0 TO NLAYERS), H(0 TO NLAYERS), D(0 TO NLAYERS), R(0 TO NLAYERS)
REDIM REFF(0 TO NLAYERS), TWT(0 TO NLAYERS), J(0 TO NLAYERS + 2)
REDIM E(0 TO NLAYERS)

FOR I = 0 TO NLAYERS - 1
PRINT "FOR LAYER # "; I; : INPUT " ENTER THE VELOCITY  "; V(I)
INPUT "                            ENTER THE THICKNESS "; H(I)
INPUT "                            ENTER THE DENSITY   "; D(I)
NEXT I
INPUT "FOR THE HALF-SPACE, ENTER THE VELOCITY "; V(NLAYERS)
INPUT "                    ENTER THE DENSITY  "; D(NLAYERS)
INPUT "                    THICKNESS (10000)  "; R$
CALL NIFDEFAULT(R$, H(NLAYERS), 10000)
ELSE
PRINT "MODEL DATA FILE NEEDS V,RHO,H (EX.  SM1.DAT)"
INPUT "ENTER THE FILE NAME OF THE MODEL DATA "; INFMOD$
INF = FREEFILE
OPEN INFMOD$ FOR INPUT AS #INF
INPUT #INF, NLAYERS
REDIM V(0 TO NLAYERS), H(0 TO NLAYERS), D(0 TO NLAYERS), R(0 TO NLAYERS)
REDIM REFF(0 TO NLAYERS), TWT(0 TO NLAYERS - 1), J(0 TO NLAYERS)
REDIM E(0 TO NP2)

FOR I = 1 TO NLAYERS
INPUT #INF, V(I), D(I), H(I)
NEXT I
CLOSE #INF
END IF

INPUT "ENTER OUTPUT FILE NAME "; OUTF$


'********************************
PRINT "NUMBER OF LAYERS "; NLAYERS
PRINT "V", "RHO", "H"
FOR I = 1 TO NLAYERS
PRINT V(I), D(I), H(I)
NEXT I

'FIND REFLECTION COEFFICIENTS
FOR N = 1 TO NLAYERS - 1
R(N) = (D(N + 1) * V(N + 1) - D(N) * V(N)) / (D(N + 1) * V(N + 1) + D(N) * V(N))
NEXT N

'FIND THE EFFECTIVE REFLECTION COEFFICIENTS
FOR N = 1 TO NLAYERS - 1
PROD = 0: TEMP = 1
FOR J = 1 TO N - 1
PROD = 1 - R(J) ^ 2
TEMP = PROD * TEMP
NEXT J
REFF(N) = TEMP * R(N)
NEXT N
PRINT "I", "R", "Reff"
FOR I = 1 TO NLAYERS - 1
PRINT I, R(I), REFF(I)
NEXT I

'FIND TWO WAY TRAVEL TIMES
COUNT = 0
FOR N = 1 TO NLAYERS - 1
SUM = 0: COUNT = COUNT + 1
FOR I = 1 TO COUNT
SUM = H(I) / V(I) * 2 + SUM
NEXT I
TWT(N) = SUM
NEXT N
PRINT "LAYER", "TWT"
FOR I = 1 TO NLAYERS - 1
PRINT I, TWT(I)
NEXT I
'INPUT DUM$
'*******************

'******************************* BP FILTER
REDIM F(0 TO 3)
'INPUT "F0 IS "; F(0)
'INPUT "F1 IS "; F(1)
'INPUT "F2 IS "; F(2)
'INPUT "F3 IS "; F(3)
DO
PRINT "ENTER THE BAND FILTER PARAMETERS"
PRINT "(1) SM1.DAT USES 60, 90, 120, 175"
PRINT "(2) SM2.DAT USES 30, 60, 100, 200"
INPUT CHOICE%
IF CHOICE% = 1 THEN
F(0) = 60: F(1) = 90: F(2) = 120: F(3) = 175         'SM1
ELSE
F(0) = 30: F(1) = 60: F(2) = 100: F(3) = 200          'SM2
END IF
LOOP WHILE CHOICE% < 1 OR CHOICE% > 2

PRINT "ENTER A PHASE (0, "; PIE$; "/2, "; PIE$; ")"
INPUT "ENTER A PHASE (0) "; R$: CALL NIFDEFAULT(R$, PHI, 0)

FMAX = 2 * F(3)

'IF FMAX < 500 THEN
'PRINT "FMAX = "; FMAX; " SO SETTING FMAX = 500 "
'FMAX = 500
'ELSE
'COLOR RED%: PRINT "WARNING !!! FMAX > 500, CAN'T SET FMAX = 500": COLOR GREEN%
'END IF

DELTAT = 1 / (2 * FMAX)
CALL LEASTDELTA(F(), LDELTAF)
PRINT "LEAST DELTAF IS "; LDELTAF; " AND 1/LDELTAF IS "; 1 / LDELTAF
FDELF = (F(1) - F(0)) / 10
PRINT "FDELF = (F(1) - F(0)) / 10 = "; FDELF; " 1/FDELF = "; 1 / FDELF
PRINT "FDELTAF = LEASTDELTAF/10 = "; LDELTAF / 10; " 1/(LDELTAF/10) = "; 1 / (LDELTAF / 10)
INPUT "ENTER THE DELTAF "; DELTAF

NP = (2 * FMAX / DELTAF)
PRINT "NP = 2*FMAX/DELTAF = "; NP: INPUT DUM$

'TMAX = CINT(TWT(NLAYERS - 1) / DELTAT + 20)           'TMAX IN MILLISECONDS
'PTMAX = TMAX
'IF TMAX > NP THEN
'NP = TMAX
'ELSE
'TMAX = NP
'END IF

CALL GETPOWEROFTWO(NP, M, NP2)
DELTAF = 1 / (NP2 * DELTAT)

A0 = -1 * LOG(.01) / ((NP2 / 4 * DELTAT) ^ 2)
TAUZERO = (NP2 / 2) * DELTAT
TIMESHIFT = NP2 / 2
START% = 0


PRINT "NP IS "; NP; " NP2 IS "; NP2
PRINT "FMAX "; FMAX; " DELTAT "; DELTAT; " 1/DELTAT = "; 1 / DELTAT
PRINT "DELTAF "; DELTAF; " 1/DELTATF = "; 1 / DELTAF

PRINT "TMAX = "; TMAX
PRINT "TAUZERO IS "; TAUZERO

UPPER = NP2

PRINT "UPPER IS "; UPPER
INPUT DUM$
REDIM TIME(START% TO 2 * UPPER)
REDIM PHASEANGLE(START% TO UPPER)
REDIM E2(START% TO UPPER), FREQ(START% TO UPPER)
REDIM RA(START% TO UPPER), IA(START% TO UPPER)
REDIM OMEGA(START% TO UPPER), FILT(START% TO UPPER)

FOR I = 0 TO NP2 - 1
PHASEANGLE(I) = PHI
NEXT I

FOR I = 0 TO UPPER - 1
TIME(I) = I * DELTAT
NEXT I
FOR I = 0 TO NP2 - 1
FREQ(I) = I * DELTAF
NEXT I
'*********************

FOR I = 0 TO NP2 - 1

IF FREQ(I) < F(0) OR FREQ(I) > F(3) THEN
FILT(I) = 0
ELSEIF FREQ(I) < F(1) THEN
FILT(I) = SIN((2 * PI * (FREQ(I) - F(0))) / (4 * (F(1) - F(0))))
ELSEIF FREQ(I) <= F(2) THEN
FILT(I) = 1
ELSE
FILT(I) = COS((2 * PI * (FREQ(I) - F(2))) / (4 * (F(3) - F(2))))
END IF
NEXT I
FPOINTS$ = STR$(F(0)) + " " + STR$(F(1)) + " " + STR$(F(2)) + " " + STR$(F(3))
FTTITLE$ = "BAND PASS FILTER" + FPOINTS$
FSTITLE$ = "BAND PASS FILTERED SIGNAL " + FPOINTS$
'*******************************

FOR I = 0 TO NP2 - 1
OMEGA(I) = FREQ(I) * 2 * PI
PHASEANGLE(I) = PHASEANGLE(I) + OMEGA(I) * TAUZERO
NEXT I

'******************** e(t) SPIKE O' GRAM

'PRINT "TWT(I) IN MILLISECONDS"
'FOR I = 1 TO NLAYERS - 1         'CONVERT SECONDS TO INTEGER MILLISECONDS
'TWT(I) = TWT(I) * 1 / DELTAT
'TWT(I) = CINT(TWT(I))
'PRINT TWT(I)
'NEXT I
'INPUT DUM$
'********************
'FIND E(J(J))

'*********************** NEW WAY

REDIM DT(0 TO NLAYERS), DN(0 TO NLAYERS), J(0 TO 1000), E(0 TO 1000)
FOR I = 1 TO NLAYERS
DT(I) = H(I) / V(I)
NEXT I

COUNT = 1
FOR I = 1 TO NLAYERS - 1
TN = H(I) / V(I) * 2
IF I = 1 THEN
DN(I) = TN
ELSE
DN(I) = DN(I - 1) + TN
END IF
J(I) = (DN(I) / DT(I)) + .5
IF COUNT = 1 THEN
MAX = J(1)
ELSEIF J(I) > J(I - 1) THEN
MAX = J(I) + 1
END IF
COUNT = COUNT + 1
NEXT I
COUNT = 1

PRINT "MAX IS "; MAX: INPUT DUM$
FOR I = 0 TO MAX
IF I = CINT(J(COUNT)) THEN
E(I) = REFF(COUNT)
COUNT = COUNT + 1
ELSE
E(I) = 0
END IF
NEXT I



'******************

'REDIM J(0 TO UPPER), E(0 TO UPPER)           'TMAX IN MILLISECONDS
'TCOUNT = 0: ECOUNT = 1
'PRINT "TMAX = "; TMAX
'PRINT "I", "J(I)", "E(I)"
'FOR I = 0 TO (UPPER - 1)               'IN INTEGER MILLISECONDS
'J(I) = I
'IF ECOUNT > NLAYERS - 1 THEN : EXIT FOR
'IF I = TWT(ECOUNT) THEN
'E(I) = REFF(ECOUNT)
'ECOUNT = ECOUNT + 1
''PRINT "I IS "; I; " E(I) IS "; E(I): INPUT DUM$
'IF I = 0 THEN
'ECOUNT = ECOUNT - 1
'END IF
'ELSE
'E(I) = 0
'END IF
''PRINT I, J(I), E(I): INPUT DUM$
'NEXT I
'ECOUNT = ECOUNT - 1
'PRINT "ECOUNT IS "; ECOUNT
'*******************************************************
'FOR I = 1 TO NLAYERS - 1   'CONVERT TWT FROM MILLISECONDS BACK TO SECONDS
'TWT(I) = TWT(I) * DELTAT
'NEXT I
'*******************************************************
'INPUT "DONE WITH E(T) "; DUM$

OPEN "E.DUM" FOR OUTPUT AS #10
'FOR I = 0 TO UPPER - 1
'PRINT #10, I, ",", E(I)
FOR I = 0 TO MAX
PRINT #10, E(I)
NEXT I
CLOSE #10

P = UPPER
CALL GETMINMAX(J(), E(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "e(t)"
YTITLE$ = "e(t)"
XTITLE$ = "t in milliseconds"
PRINT "J", "E"

'XMAX = PTMAX
INPUT DUM$
CALL GRAPH(J(), E(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)
'********************
OUTF = FREEFILE
OPEN OUTF$ FOR OUTPUT AS #OUTF
PRINT "V", "H", "RHO", "R", "Reff", "TWT", "J", "E(J)"
PRINT #OUTF, "V", "H", "RHO", "R", "Reff", "TWT", "J", "E(J)"
FOR I = 1 TO NLAYERS - 1
PRINT #OUTF, V(I), ",", H(I), ",", D(I), ",", R(I), ",", REFF(I), ",", TWT(I), ",", J(I), ",", E(I)
PRINT V(I), H(I), D(I), R(I), REFF(I), TWT(I), J(I), E(I)
NEXT I
PRINT #OUTF, V(NLAYERS), ",", H(NLAYERS), ",", D(NLAYERS)
PRINT V(NLAYERS), H(NLAYERS), D(NLAYERS)

'******************************

PRINT "READY TO GRAPH AMPLITUDE  SPECTRUM (THE FILTER)"
P = NP2
CALL GETMINMAX(FREQ(), FILT(), P, XMIN, XMAX, YMIN, YMAX)
T$ = FSTITLE$ + " AMPLITUDE (FILTER)"
XMAX = FMAX
CALL GRAPH(FREQ(), FILT(), P, XMIN, XMAX, YMIN, YMAX, 0, T$, "FREQUENCY (w) Hz", "S(w)")
'************************

FOR I = 0 TO NP2 - 1  'CONVERT BACK TO RECTANGULAR BEFORE FFT
RA(I) = FILT(I) * COS(PHASEANGLE(I))
IA(I) = FILT(I) * (SIN(PHASEANGLE(I))) * -1
NEXT I

CALL FFT(NP2, M, RA(), IA())    'SENDING F DOMAIN, WILL RETURN T DOMAIN

REDIM BELL(0 TO NP2)
FOR I = 0 TO NP2 - 1
'RA(I) = RA(I) * EXP(-1 * A0 * (TIME(I) - TAUZERO) ^ 2)
TEMP = I * DELTAT
BELL(I) = EXP(-A0 * (TEMP - TAUZERO) ^ 2)
RA(I) = RA(I) * BELL(I)
NEXT I

P = NP2 - 1
PRINT "GRAPHING BELL CURVE"
CALL GETMINMAX(TIME(), BELL(), P, XMIN, XMAX, YMIN, YMAX)
FSTITLE$ = "BELL CURVE"
CALL GRAPH(TIME(), BELL(), P, XMIN, XMAX, YMIN, YMAX, 0, FSTITLE$, "TIME (SEC)", "S(t)")

CALL GETMINMAX(TIME(), RA(), P, XMIN, XMAX, YMIN, YMAX)
FSTITLE$ = "THE NEW PULSE"
CALL GRAPH(TIME(), RA(), P, XMIN, XMAX, YMIN, YMAX, 0, FSTITLE$, "TIME (SEC)", "S(t)")

OPEN "RA.DUM" FOR OUTPUT AS #30
PRINT #30, "THE PULSE"
PRINT #30, "I", "TIME", "RA(I)"
FOR I = 0 TO NP2 - 1
PRINT #30, I, ",", TIME(I), ",", RA(I)
NEXT I
CLOSE #30

'***********************
REDIM CON(0 TO UPPER * 2), RAFOLD(0 TO UPPER)      'E2(0 TO NP2)
REDIM E2(0 TO UPPER), RA2(0 TO UPPER), CONCOUNT(0 TO UPPER * 2)
FOR I = 0 TO UPPER - 1
E2(I) = E(I)
NEXT I
FOR I = 0 TO NP2 - 1
RA2(I) = RA(I)
NEXT I

CALL FOLDANDSHIFT(NP2, RA(), RAFOLD())

CT$ = "e(t) * s(t)"

CALL NEWCONV(TMAX, NP2, E2(), RA2(), RAFOLD(), CON(), CONCOUNT())
'    NEWCONV (ISNP, IGNP, S(), G(), GF(), CON(), CONCOUNT())

'*****************************************************************
PRINT "BACK IN SEISMOGRAM": INPUT DUM$

FOR I = 0 TO 2 * UPPER - 1
CONCOUNT(I) = CONCOUNT(I) - TAUZERO * 1000
NEXT I

'*****************************************
PRINT "READY TO GRAPH S(t)"
P = UPPER
CALL GETMINMAX(J(), E2(), P, XMIN, XMAX, YMIN, YMAX)
XMAX = PTMAX
TITLE$ = "e(t)"
XTITLE$ = "time (milliseconds)"
CALL GRAPH(J(), E2(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE PULSE"
P = UPPER
CALL GETMINMAX(TIME(), RA2(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "THE PULSE"
'YTITLE$ = "g(t)"
XTITLE$ = "time (milliseconds)"
CALL GRAPH(TIME(), RA2(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE FOLDED SIGNAL"
P = UPPER
CALL GETMINMAX(TIME(), RAFOLD(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "FOLDED PULSE"
'YTITLE$ = "G(t) FOLDED"
XTITLE$ = "time (milliseconds)"
'XMAX = PTMAX
CALL GRAPH(TIME(), RAFOLD(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)


PRINT "READY TO GRAPH THE CONVOLUTION"
'LINE INPUT "ENTER THE TITLE "; TITLE$
TITLE$ = "p(t) = e(t) * s(t)"
P = UBOUND(CONCOUNT) - LBOUND(CONCOUNT)
CALL GETMINMAX(CONCOUNT(), CON(), P, XMIN, XMAX, YMIN, YMAX)
XTITLE$ = "time (milliseconds)"
XMAX = (UBOUND(CONCOUNT) - LBOUND(CONCOUNT)) / 2
XMIN = 0
CALL GRAPH(CONCOUNT(), CON(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)


'********************************************************
PRINT "WRITING OUTPUT TO "; OUTF$
PRINT #OUTF, "I", "E(T)", "PULSE()"
FOR I = 0 TO NP2
PRINT #OUTF, I, ",", E2(I), ",", RA(I)
NEXT I
PRINT #OUTF, "I", "CON()"
FOR I = LBOUND(CON) TO UBOUND(CON) - 1
PRINT #OUTF, I, ",", CON(I)
NEXT I
CLOSE #OUTF
ERASE CON, E2, RA, IA, RAFOLD, PHASEANGLE, TIME, OMEGA, FREQ, FILT
ERASE BELL, V, H, D, R, REFF, TWT, J, E, CONCOUNT, RA2


END SUB

SUB SEISMOGRAM (NEWMODEL)
'***********************  GET E(t)
DEL$ = CHR$(127)
PIE$ = CHR$(227)
PHI$ = CHR$(237)
RHO$ = CHR$(235)

NLAYERS = 0:

PRINT "PHI "; PHI$; " DEL "; DEL$; " PI "; PIE$; " RHO "; RHO$

IF NEWMODEL THEN
'INPUT "ENTER DELTA T IN SECONDS "; DELTAT
INPUT "ENTER THE NUMBER OF LAYERS INCLUDING THE HALF SPACE "; NLAYERS
REDIM V(0 TO NLAYERS), H(0 TO NLAYERS), D(0 TO NLAYERS), R(0 TO NLAYERS)
REDIM REFF(0 TO NLAYERS), TWT(0 TO NLAYERS), J(0 TO NLAYERS + 2)
REDIM E(0 TO NLAYERS)

FOR I = 0 TO NLAYERS - 1
PRINT "FOR LAYER # "; I; : INPUT " ENTER THE VELOCITY  "; V(I)
INPUT "                            ENTER THE THICKNESS "; H(I)
INPUT "                            ENTER THE DENSITY   "; D(I)
NEXT I
INPUT "FOR THE HALF-SPACE, ENTER THE VELOCITY "; V(NLAYERS)
INPUT "                    ENTER THE DENSITY  "; D(NLAYERS)
INPUT "                    THICKNESS (10000)  "; R$
CALL NIFDEFAULT(R$, H(NLAYERS), 10000)
ELSE
PRINT "MODEL DATA FILE NEEDS V,RHO,H (EX.  SM1.DAT)"
FILES "*.DAT"
INPUT "ENTER THE FILE NAME OF THE MODEL DATA "; INFMOD$
INF = FREEFILE
OPEN INFMOD$ FOR INPUT AS #INF
INPUT #INF, NLAYERS
REDIM V(0 TO NLAYERS), H(0 TO NLAYERS), D(0 TO NLAYERS), R(0 TO NLAYERS)
REDIM REFF(0 TO NLAYERS), TWT(0 TO NLAYERS - 1), J(0 TO NLAYERS)
REDIM E(0 TO NP2)

FOR I = 1 TO NLAYERS
INPUT #INF, V(I), D(I), H(I)
NEXT I
CLOSE #INF
END IF

INPUT "ENTER OUTPUT FILE NAME "; OUTF$

INPUT "ENTER ADDITIONAL TITLE "; COMMONTITLE$
'********************************
PRINT "NUMBER OF LAYERS "; NLAYERS
PRINT "V", "RHO", "H"
FOR I = 1 TO NLAYERS
PRINT V(I), D(I), H(I)
NEXT I

'FIND REFLECTION COEFFICIENTS
FOR N = 1 TO NLAYERS - 1
R(N) = (D(N + 1) * V(N + 1) - D(N) * V(N)) / (D(N + 1) * V(N + 1) + D(N) * V(N))
NEXT N

'FIND THE EFFECTIVE REFLECTION COEFFICIENTS
FOR N = 1 TO NLAYERS - 1
PROD = 0: TEMP = 1
FOR J = 1 TO N - 1
PROD = 1 - R(J) ^ 2
TEMP = PROD * TEMP
NEXT J
REFF(N) = TEMP * R(N)
NEXT N
PRINT "I", "R", "Reff"
FOR I = 1 TO NLAYERS - 1
PRINT I, R(I), REFF(I)
NEXT I

'FIND TWO WAY TRAVEL TIMES
COUNT = 0
FOR N = 1 TO NLAYERS - 1
SUM = 0: COUNT = COUNT + 1
FOR I = 1 TO COUNT
SUM = H(I) / V(I) * 2 + SUM
NEXT I
TWT(N) = SUM
NEXT N
PRINT "LAYER", "TWT"
FOR I = 1 TO NLAYERS - 1
PRINT I, TWT(I)
NEXT I
'INPUT DUM$
'*******************     MAKING LOG FILE
INPUT "APPEND TO LOG FILE SMOD.LOG (Y) "; YN$
IF YN$ = "" THEN
YN$ = "Y"
END IF
IF YN$ = "Y" THEN
OL = FREEFILE
OPEN "SMOD.LOG" FOR APPEND AS #OL
PRINT #OL, INFMOD$
PRINT #OL, "V", "RHO", "H"
FOR I = 1 TO NLAYERS
PRINT #OL, V(I), ",", D(I), ",", H(I)
NEXT I

PRINT #OL, "R", "Reff", "TWT(sec)", "TWT(msec)"
FOR N = 1 TO NLAYERS - 1
PRINT #OL, R(N), ",", REFF(N), ",", TWT(N), ",", TWT(N) * 1000
NEXT N
PRINT #OL, ""
CLOSE #OL
END IF

'******************************* BP FILTER
REDIM F(0 TO 3)
DO
PRINT "ENTER THE BAND FILTER PARAMETERS"
PRINT "(1) SM1.DAT USES 60, 90, 120, 175"
PRINT "(2) SM2.DAT USES 30, 60, 100, 200"
PRINT "(3) SM2.DAT USEING 20, 40, 50, 100"
PRINT
PRINT "(4) INPUT BAND PASS FILTER "
INPUT CHOICE%
LOOP WHILE CHOICE% < 1 OR CHOICE% > 4
SELECT CASE CHOICE%
CASE 1:
F(0) = 60: F(1) = 90: F(2) = 120: F(3) = 175         'SM1
CASE 2:
F(0) = 30: F(1) = 60: F(2) = 100: F(3) = 200          'SM2
CASE 3:
F(0) = 20: F(1) = 40: F(2) = 50: F(3) = 100
CASE 4:
INPUT "F0 IS "; F(0)
INPUT "F1 IS "; F(1)
INPUT "F2 IS "; F(2)
INPUT "F3 IS "; F(3)
END SELECT


PRINT "ENTER A PHASE (0, "; PIE$; "/2, "; PIE$; ")"
PRINT "ENTER A PHASE ("; PIE$; "/2) ": INPUT R$: CALL NIFDEFAULT(R$, PHI, 1.57)
PRINT "PHASE IS "; PHI
FMAX = 2 * F(3)
IF FMAX < 500 THEN
PRINT "FMAX = "; FMAX; " SO SETTING FMAX = 500 "
FMAX = 500
ELSE
COLOR RED%: PRINT "WARNING !!! FMAX > 500, CAN'T SET FMAX = 500"
PRINT "FMAX WILL BE SET TO "; FMAX
END IF

DELTAT = 1 / (2 * FMAX)
CALL LEASTDELTA(F(), DELTAF)
PRINT "LEAST DELTAF IS "; DELTAF; " AND 1/DELTAF IS "; 1 / DELTAF
NP = (2 * FMAX / DELTAF)
TMAX = CINT(TWT(NLAYERS - 1) / DELTAT + 20)           'TMAX IN MILLISECONDS
PTMAX = TMAX
IF TMAX > NP THEN
NP = TMAX
ELSE
TMAX = NP
END IF

CALL GETPOWEROFTWO(NP, M, NP2)
DELTAF = 1 / (NP2 * DELTAT)

A0 = -1 * LOG(.01) / ((NP2 / 4 * DELTAT) ^ 2)
'TAUZERO = (NP2 / 2) * DELTAT
'TAUZERO = NP * DELTAT
TAUZERO = (NP2 - 1 - NP) * DELTAT
TIMESHIFT = NP2 / 2
START% = 0


PRINT "NP IS "; NP; " NP2 IS "; NP2
PRINT "FMAX "; FMAX; " DELTAT "; DELTAT; " 1/DELTAT = "; 1 / DELTAT
PRINT "DELTAF "; DELTAF; " 1/DELTATF = "; 1 / DELTAF

PRINT "TMAX = "; TMAX
PRINT "TAUZERO IS "; TAUZERO

UPPER = NP2
PRINT "UPPER IS "; UPPER
INPUT DUM$
REDIM TIME(START% TO 2 * UPPER)
REDIM PHASEANGLE(START% TO UPPER)
REDIM E2(START% TO UPPER), FREQ(START% TO UPPER)
REDIM RA(START% TO UPPER), IA(START% TO UPPER)
REDIM OMEGA(START% TO UPPER), FILT(START% TO UPPER)

FOR I = 0 TO NP2 - 1
PHASEANGLE(I) = PHI
NEXT I

FOR I = 0 TO UPPER - 1
TIME(I) = I * DELTAT
NEXT I
FOR I = 0 TO NP2 - 1
FREQ(I) = I * DELTAF
NEXT I
'*********************

FOR I = 0 TO NP2 - 1

IF FREQ(I) < F(0) OR FREQ(I) > F(3) THEN
FILT(I) = 0
ELSEIF FREQ(I) < F(1) THEN
FILT(I) = SIN((2 * PI * (FREQ(I) - F(0))) / (4 * (F(1) - F(0))))
ELSEIF FREQ(I) <= F(2) THEN
FILT(I) = 1
ELSE
FILT(I) = COS((2 * PI * (FREQ(I) - F(2))) / (4 * (F(3) - F(2))))
END IF
NEXT I
FPOINTS$ = STR$(F(0)) + " " + STR$(F(1)) + " " + STR$(F(2)) + " " + STR$(F(3))
FTTITLE$ = "BAND PASS FILTER" + FPOINTS$
FSTITLE$ = "BAND PASS FILTERED SIGNAL " + FPOINTS$
'*********************************
FOR I = 0 TO NP2 - 1
OMEGA(I) = FREQ(I) * 2 * PI
PHASEANGLE(I) = PHASEANGLE(I) + OMEGA(I) * TAUZERO
NEXT I

'******************** e(t) SPIKE O' GRAM

PRINT "TWT(I) IN MILLISECONDS"
FOR I = 1 TO NLAYERS - 1         'CONVERT SECONDS TO INTEGER MILLISECONDS
TWT(I) = TWT(I) * 1 / DELTAT
TWT(I) = CINT(TWT(I))
PRINT TWT(I)
NEXT I
INPUT DUM$
'********************
'FIND E(J(J))

REDIM TSCALE(0 TO UPPER), E(0 TO UPPER)       'TMAX IN MILLISECONDS
TCOUNT = 0: ECOUNT = 1
PRINT "TMAX = "; TMAX
PRINT "I", "TSCALE(I)", "E(I)"
FOR I = 0 TO (UPPER - 1)             'IN INTEGER MILLISECONDS
TSCALE(I) = I
IF ECOUNT > NLAYERS - 1 THEN : EXIT FOR
IF I = TWT(ECOUNT) THEN
E(I) = REFF(ECOUNT)
ECOUNT = ECOUNT + 1
'PRINT "I IS "; I; " E(I) IS "; E(I): INPUT DUM$
IF I = 0 THEN
ECOUNT = ECOUNT - 1
END IF
ELSE
E(I) = 0
END IF
'PRINT I, TSCALE(I), E(I): INPUT DUM$
NEXT I
ECOUNT = ECOUNT - 1
PRINT "ECOUNT IS "; ECOUNT
'*******************************************************
FOR I = 1 TO NLAYERS - 1   'CONVERT TWT FROM MILLISECONDS BACK TO SECONDS
TWT(I) = TWT(I) * DELTAT
NEXT I
'*******************************************************
'INPUT "DONE WITH E(T) "; DUM$

P = UPPER
CALL GETMINMAX(TSCALE(), E(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "e(t)"
YTITLE$ = "e(t)"
XTITLE$ = "t in milliseconds"
PRINT "J", "E"
XMAX = PTMAX
CALL GRAPH(TSCALE(), E(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)
'********************
OUTF = FREEFILE
OPEN OUTF$ FOR OUTPUT AS #OUTF
PRINT "V", "H", "RHO", "R", "Reff", "TWT", "J", "E(J)"
PRINT #OUTF, "V", "H", "RHO", "R", "Reff", "TWT", "J", "E(J)"
FOR I = 1 TO NLAYERS - 1
PRINT #OUTF, V(I), ",", H(I), ",", D(I), ",", R(I), ",", REFF(I), ",", TWT(I), ",", J(I), ",", E(I)
PRINT V(I), H(I), D(I), R(I), REFF(I), TWT(I), J(I), E(I)
NEXT I
PRINT #OUTF, V(NLAYERS), ",", H(NLAYERS), ",", D(NLAYERS)
PRINT V(NLAYERS), H(NLAYERS), D(NLAYERS)

'******************************

PRINT "READY TO GRAPH AMPLITUDE  SPECTRUM (THE FILTER)"
P = NP2
CALL GETMINMAX(FREQ(), FILT(), P, XMIN, XMAX, YMIN, YMAX)
T$ = FSTITLE$ + " AMPLITUDE (FILTER)"
XMAX = FMAX
CALL GRAPH(FREQ(), FILT(), P, XMIN, XMAX, YMIN, YMAX, 0, T$, "FREQUENCY (w) Hz", "S(w)")
'************************

FOR I = 0 TO NP2 - 1  'CONVERT BACK TO RECTANGULAR BEFORE FFT
RA(I) = FILT(I) * COS(PHASEANGLE(I))
IA(I) = FILT(I) * (SIN(PHASEANGLE(I))) * -1
NEXT I

CALL FFT(NP2, M, RA(), IA())    'SENDING F DOMAIN, WILL RETURN T DOMAIN

REDIM BELL(0 TO NP2)
FOR I = 0 TO NP2 - 1
'RA(I) = RA(I) * EXP(-1 * A0 * (TIME(I) - TAUZERO) ^ 2)
TEMP = I * DELTAT
BELL(I) = EXP(-A0 * (TEMP - TAUZERO) ^ 2)
RA(I) = RA(I) * BELL(I)
NEXT I

P = NP2 - 1
PRINT "GRAPHING BELL CURVE"
CALL GETMINMAX(TIME(), BELL(), P, XMIN, XMAX, YMIN, YMAX)
FSTITLE$ = "BELL CURVE"
CALL GRAPH(TIME(), BELL(), P, XMIN, XMAX, YMIN, YMAX, 0, FSTITLE$, "TIME (SEC)", "S(t)")

CALL GETMINMAX(TIME(), RA(), P, XMIN, XMAX, YMIN, YMAX)
FSTITLE$ = "THE NEW PULSE"
CALL GRAPH(TIME(), RA(), P, XMIN, XMAX, YMIN, YMAX, 0, FSTITLE$, "TIME (SEC)", "S(t)")

'***********************
REDIM CON(0 TO UPPER * 2), RAFOLD(0 TO UPPER)      'E2(0 TO NP2)
REDIM E2(0 TO UPPER), RA2(0 TO UPPER), CONCOUNT(0 TO UPPER * 2)
FOR I = 0 TO UPPER - 1
E2(I) = E(I)
NEXT I
FOR I = 0 TO NP2 - 1
RA2(I) = RA(I)
NEXT I

CALL FOLDANDSHIFT(NP2, RA(), RAFOLD())
CT$ = "e(t) * s(t)"

CALL NEWCONV(TMAX, NP2, E2(), RA2(), RAFOLD(), CON(), CONCOUNT())
'    NEWCONV (ISNP, IGNP, S(), G(), GF(), CON(), CONCOUNT())

'*****************************************************************
PRINT "TAUZERO IS "; TAUZERO
PRINT "BACK IN SEISMOGRAM": INPUT DUM$
TAUZERO = TAUZERO - DELTAT
FOR I = 0 TO 2 * UPPER - 1
CONCOUNT(I) = CINT(CONCOUNT(I) - TAUZERO * 1 / DELTAT)
NEXT I
TAUZERO = CINT(TAUZERO * 1 / DELTAT)
'*****************************************
PRINT "READY TO GRAPH S(t)"
P = UPPER
CALL GETMINMAX(TSCALE(), E2(), P, XMIN, XMAX, YMIN, YMAX)
XMAX = PTMAX
TITLE$ = "e(t)" + FPOINTS$
XTITLE$ = "time (milliseconds)"
CALL GRAPH(TSCALE(), E2(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE PULSE"
P = UPPER
CALL GETMINMAX(TIME(), RA2(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "THE PULSE" + FPOINTS$
'YTITLE$ = "g(t)"
XTITLE$ = "time (milliseconds)"
CALL GRAPH(TIME(), RA2(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)

PRINT "READY TO GRAPH THE FOLDED SIGNAL"
P = UPPER
CALL GETMINMAX(TIME(), RAFOLD(), P, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "FOLDED PULSE" + FPOINTS$
'YTITLE$ = "G(t) FOLDED"
XTITLE$ = "time (milliseconds)"
'XMAX = PTMAX
CALL GRAPH(TIME(), RAFOLD(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)


PRINT "READY TO GRAPH THE CONVOLUTION"
'LINE INPUT "ENTER THE TITLE "; TITLE$
TITLE$ = "p(t) = e(t) * s(t) " + FPOINTS$
P = UBOUND(CONCOUNT) - LBOUND(CONCOUNT)
CALL GETMINMAX(CONCOUNT(), CON(), P, XMIN, XMAX, YMIN, YMAX)
XTITLE$ = "time (milliseconds)"
XMAX = (UBOUND(CONCOUNT) - LBOUND(CONCOUNT)) / 2
PRINT "XMAX "; XMAX
'XMIN = 0
XMIN = CONCOUNT(TAUZERO)
CALL GRAPH(CONCOUNT(), CON(), P, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)


'********************************************************
PRINT "WRITING OUTPUT TO "; OUTF$
PRINT #OUTF, "I", "E(T)", "PULSE()"
FOR I = 0 TO NP2
PRINT #OUTF, I, ",", E2(I), ",", RA(I)
NEXT I
PRINT #OUTF, "I", "CONCOUNT", "CON()"
FOR I = LBOUND(CON) TO UBOUND(CON) - 1
PRINT #OUTF, I, ",", CONCOUNT(I), ",", CON(I)
NEXT I
CLOSE #OUTF
ERASE CON, E2, RA, IA, RAFOLD, PHASEANGLE, TIME, OMEGA, FREQ, FILT
ERASE BELL, V, H, D, R, REFF, TWT, J, E, CONCOUNT, RA2

END SUB

SUB SEISMOGRAMMENU
'ARTIFICIAL SEISMOGRAM
DO
COLOR 2
        'DISPLAY MENU
CLS 0
PRINT "ARTIFICIAL SEISMOGRAM"
PRINT ""
PRINT "         (1) ENTER A NEW MODEL"
PRINT "         (2) INPUT A PREVIOUS MODEL"
PRINT "         (3) "
PRINT "         (4)"
PRINT ""
PRINT "         (5) EXIT TO DOS SHELL (TYPE 'EXIT' TO RETURN TO PROGRAM)"
PRINT ""
PRINT "         (6) CONVOLUTION"
PRINT
PRINT "         (7)"
PRINT ""
PRINT "         (8) VER 2. ENTER A NEW MODEL"
PRINT ""
PRINT "         (9) VER 2. INPUT A MODEL"
PRINT
PRINT "        (10) QUIT THIS MENU"
INPUT "ENTER A NUMBER TO INDICATE YOUR CHOICE "; CHOICE%
PRINT ""
'RESPOND TO CHOICE
SELECT CASE CHOICE%
        CASE 1:
        'INPUT A NEW MODEL
            NEWMODEL = TRUE
            CALL SEISMOGRAM(NEWMODEL)
        CASE 2:
        'INPUT AN EXISTING MODEL
        CALL SEISMOGRAM(NEWMODEL)

        CASE 3:

        CASE 4:




        CASE 5: COLOR 9
        PRINT "TYPE EXIT TO RETURN TO PROGRAM"
                SHELL
                COLOR 2

        CASE 6:  CALL NEWCONVOLUTION


        CASE 7:

        CASE 8:
        'INPUT A NEW MODEL
            NEWMODEL = TRUE
            CALL SEISMO2(NEWMODEL)



        CASE 9:
        'INPUT AN EXISTING MODEL
        CALL SEISMO2(NEWMODEL)


        CASE 10: DONE = TRUE


END SELECT       'SELECT CASE CHOICE
LOOP UNTIL DONE




END SUB

SUB SINCLPF (NP, DELT, F())
PRINT "LOW PASS FILTER BY CONVOLUTION"
INPUT "ENTER F0 (40) "; R$: CALL NIFDEFAULT(R$, F0, 40)
OUF = FREEFILE
OPEN "F40SINC.DAT" FOR OUTPUT AS #OUF

REDIM FCOUNT(0 TO 2 * NP)
F(0) = 0
K = 0
DT = 0
FOR I = DELT * -(NP / 2) TO DELT * ((NP / 2)) STEP DELT
DT = I
IF DT <> 0 THEN
F(K) = 2 * F0 * SIN(2 * PI * F0 * DT) / (2 * PI * F0 * DT)
ELSE
F(K) = 2 * F0
END IF
PRINT #OUF, F(K)
FCOUNT(K) = K
K = K + 1
NEXT I
CLOSE #OUF
FTOTAL = K
CALL GETMINMAX(FCOUNT(), F(), FTOTAL, XMIN, XMAX, YMIN, YMAX)
TITLE$ = "f(t)"
YTITLE$ = "f(t)"
XTITLE$ = "t  DELTAT = " + STR$(DELT)
CALL GRAPH(FCOUNT(), F(), FTOTAL, XMIN, XMAX, YMIN, YMAX, 0, TITLE$, XTITLE$, YTITLE$)
ERASE FCOUNT
END SUB

SUB WAITFORKEY
DO: LOOP UNTIL INKEY$ <> ""
END SUB